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Abstract 

The  influence  of  particle-size  distribution  on  the  stress  transfer  mecha¬ 
nisms  occurring  in  bimodal  granular  materials  was  investigated  in  this 
study.  Efficient  stress  transfer  is  achieved  when  all  particles  within  an 
assembly  contribute  to  the  resistance  of  applied  loads.  Both  physical  and 
numerical  testing  methods  were  employed  including  consolidated- 
undrained  triaxial  testing  and  numerical  modeling  using  the  discrete 
element  method  (DEM).  Materials  investigated  included  ideal  spheres  of 
stainless  steel,  polypropylene,  and  natural  materials  consisting  of  sand  and 
silty-clay.  Identification  of  critical  mixture  proportions  according  to 
percolation  theory  was  attempted  by  examining  the  macro-  and  micro¬ 
mechanical  response  of  bimodal  distributions  that  were  dominated  by  one 
fraction  or  were  approaching  the  percolation  threshold.  Research  results 
were  inconclusive  in  validating  the  original  hypothesis  when  comparing 
the  numerical  results  to  laboratory  experiments.  However,  a  relationship 
between  coordination  number  and  assembly  stiffness  was  observed. 
Factors  contributing  to  increased  connectivity  can  also  contribute  to 
increased  stiffness  and  shear  strength. 
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1  Introduction 

The  design  and  evaluation  of  engineered  granular  media  requires  an 
understanding  of  the  magnitude  and  distribution  of  stresses  induced  from 
application  specific  loading.  Granular  or  particulate  media  includes  soils, 
powders,  agricultural  grains,  and  industrial  materials.  The  stress  state  of  a 
granular  assembly  can  be  used  to  explain  and  predict  engineering 
behavior.  The  response  of  particulate  systems  to  externally  applied  loads  is 
typically  explained  using  constitutive  theories,  such  as  the  critical  state 
theory  for  soils,  based  on  laboratory  studies  of  macro-scale  behavior. 

These  methods  rely  upon  measurements  of  boundary  conditions  to 
estimate  the  stresses  within  the  specimen  due  to  the  difficulty  of  collecting 
internal  measurements.  Constitutive  methods  also  require  treating  the 
particulate  system  as  homogeneous,  elastic,  and  isotropic.  Most 
particulate  systems,  in  particular  soils,  deviate  from  these  idealizations 
often  resulting  in  stress  calculation  errors  ranging  from  25  to  30  percent 
(Das  2010).  Accurate  prediction  of  the  mechanical  response  of  particulate 
systems  requires  an  understanding  of  the  stress  transmission  mechanisms 
within  an  assembly  of  particles  without  the  constraints  imposed  by 
idealized  constitutive  theory  based  methods. 

A  number  of  researchers  have  begun  examining  the  engineering  response  of 
particle  systems  using  discrete  numerical  methods.  Numerical  models 
provide  data  on  the  particle-level  mechanics,  such  as  displacement,  rota¬ 
tion,  and  inter-  and  intra-particle  forces,  which  cannot  be  obtained  through 
traditional  laboratory  instrumentation.  Finite,  particle  based  methods  yield 
results  more  consistent  with  laboratory  observations  as  compared  to 
continuum  methods  (Armstrong  and  Dunlap  1966).  Mathematical  models 
such  as  the  discrete  element  method  (DEM)  allow  for  simulation  of 
traditional  laboratory  testing  methods  while  also  providing  data  on  the 
micro-  or  particle-level  mechanics  driving  the  system.  The  DEM  is  an 
explicit  numerical  approach  in  which  the  interaction  of  particles  is 
monitored  contact  by  contact  and  the  motion  is  monitored  particle  by 
particle  (Cundall  and  Strack  1979).  Data  from  the  DEM  can  be  used  to 
investigate  the  micro-level  physics,  which  control  the  observable,  macro¬ 
mechanical  performance  of  the  material.  This  research  seeks  to  use  a  DEM 
model  to  investigate  stress  transfer  mechanisms  by  quantifying  the  influ¬ 
ence  of  particle-size  distribution  on  engineering  performance.  The  principal 
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advantage  of  using  numerical  simulations  is  that  all  inter-particle  data  can 
be  extracted  from  the  assembly  under  investigation.  Additionally,  the 
influence  of  micromechanical  properties  can  be  more  readily  determined  as 
compared  to  physical  experiments  (Bathurst  and  Rothenburg  1988). 

1.1  Research  objective 

The  transmission  of  stresses  between  particles  is  based  upon  a  complex 
relationship  determined  by  the  particle-size  distribution,  loading  history, 
fabric  (assembly  configuration),  and  particle  characteristics  (stiffness, 
shape,  roughness  etc.).  Efficient  stress  transfer  occurs  when  all  of  the 
particles  in  a  consolidated  assembly  participate  in  resisting  applied  loads. 
The  objective  of  this  research  is  to  investigate  the  influence  of  particle-size 
distribution  on  the  stress  transfer  mechanisms  occurring  in  binary  particle 
systems.  Many  gap-graded  materials  commonly  used  in  engineering  can 
be  approximated  using  binary  systems  (Peters  and  Berney  2010).  This 
effort  seeks  to  describe  the  engineering  behavior  of  binary  blends  of  parti¬ 
cles  by  examining  the  micro-  and  macro-mechanical  response  of  the 
systems  to  loading  using  both  physical  and  numerical  experiments. 

1.2  Research  approach 

The  objective  of  this  research  will  be  accomplished  by  examining  the 
micro-mechanical  data  obtained  from  a  DEM  model  to  quantify  the 
influence  of  particle-size  distribution  on  macro-scale  engineering 
performance.  The  DEM  model  developed  by  Cundall  and  Strack  (1979)  as 
implemented  by  Horner  et.al  (2001)  will  be  used  to  simulate  the  triaxial 
compression  testing  of  consolidated  cubical  specimens  of  particles. 
Systems  of  consolidated  spheres  in  rhombic  packing  configurations  will  be 
generated  using  identical  methods  to  eliminate  the  influence  of  other 
factors  contributing  to  stress  transfer  mechanisms.  The  DEM  model  will 
provide  data  on  the  inter-  and  intra-particle  forces,  allowing  for 
examination  of  the  distribution  of  stresses  and  connectivity  evolution 
amongst  individual  particles  and  between  the  discrete  size  fractions. 

Additionally,  the  macro-mechanical  properties  of  representative  physical 
specimens  will  be  determined  by  laboratory  triaxial  compression  testing 
conducted  in  accordance  with  ASTM  International  standard  test  method 
D4767,  Standard  test  method  for  consolidated  undrained  triaxial 
compression  test  for  cohesive  soils  (ASTM  2004).  The  triaxial  test  is  the 
most  commonly  used  test  for  determining  the  constitutive  response  of 
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geotechnical  materials.  Although  recent  technical  advancements  in  data 
collection  have  allowed  measurement  of  some  internal  response  within  a 
triaxial  test  specimen  including  the  evolution  of  localizations,  some  data 
are  unattainable  through  physical  testing  including  inter-particle  forces 
(Cui  et.  al  2007). 

The  distribution  of  particles  within  discrete  size  fractions  will  be  varied  to 
define  the  influence  of  gradation  on  macro-mechanical  properties.  Both 
ideal  spherical  and  naturally  occurring  particles  will  be  examined  in  this 
investigation.  Ideal  particles  will  consist  of  discrete  binary  blends  of  stain¬ 
less  steel  and  polypropylene  spheres.  The  particle-size  distribution  of 
many  important  engineering  applications  of  granular  materials  can  be 
approximated  as  binary  mixtures  including  asphalt  pavement,  retaining 
structures,  and  dam  and  levee  layers  (Peters  and  Berney  2010).  Idealizing 
soil  grains  as  spherical  particles  departs  from  reality  however;  the 
idealization  is  sustainable  when  it  provides  useful  information  about 
material  response  (Bardet  1998).  Naturally  occurring  particles  will  be 
investigated  using  data  provided  in  Peters  and  Berney  (2010).  The  authors 
used  bi-modal  grain-size  distributions  of  sand  and  silty-clay  particles  to 
examine  the  theory  of  percolation.  Future  research  is  intended  to  consider 
real  world  particles  represented  by  sand  and  silty-clay  grain  fractions. 

1.3  Organization  of  report 

Chapter  1  presents  an  introduction  to  the  research  topic.  Chapter  2  of  this 
document  provides  essential  supporting  background  information  used  in 
development  of  the  problem  statement.  The  DEM  numerical  model  is 
presented  in  Chapter  3  including  the  key  features,  history  of  development, 
and  calibration  parameters.  Chapter  4  presents  the  details  of  the 
numerical  simulations  run  as  part  of  this  investigation,  and  Chapter  5 
presents  the  details  of  the  laboratory  testing.  Chapter  6  presents  analysis 
of  results  and  relevant  discussion.  A  summary  of  the  conducted  research 
and  key  conclusions  resulting  from  the  investigation  is  presented  in 
Chapter  7.  Finally,  a  list  of  referenced  documents  is  presented. 
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2  Background 

Any  investigation  of  the  behavior  of  granular  materials  in  engineering 
applications  requires  knowledge  of  the  stress  and  strain  distributions 
occurring  in  the  material  (Oda  and  Iwashita  1999).  Traditional  methods 
for  determining  stresses  (i.e.,  Boussinesq’s  equation  for  soils)  require 
treating  the  particulate  system  as  an  idealized  homogeneous,  elastic,  and 
isotropic  system  (Das  2010).  Most  particulate  systems,  in  particular  soils, 
deviate  from  these  idealizations  often  resulting  in  calculation  errors  of 
25  to  30  percent  (Das  2010).  The  assumed  stresses  in  a  particulate  system 
under  a  given  loading  configuration  are  the  primary  drivers  in  the  design 
and  evaluation  of  engineering  systems. 

Micromechanics  explains  the  relationship  between  externally  applied 
stresses  and  strains  and  internal  forces  and  displacements  (Sitharam  et  al. 
2002).  The  micromechanics  of  granular  media  is  a  topic  of  interest  for  a 
number  of  technical  fields  including  civil  engineering,  chemistry, 
metallurgy,  pharmacy,  and  electrical  engineering.  Relationships  between 
stress  and  strain  can  be  examined  phenomenologically  or  structurally.  The 
traditional  phenomenological  approach  relates  measured  mechanical 
properties  to  observed  behavior  and  is  analogous  with  constitutive  theory 
based  methods.  A  structural  approach  examines  the  mechanical  behavior 
based  upon  the  interactions  of  fundamental  constitutive  units.  Finite 
methods  are  representative  of  the  structural  approach.  Feda  (1982) 
comments  that  the  structural  approach  is  required  to  understand  the 
physical  fundamentals  of  the  mechanical  behavior  of  particulate  materials. 

2.1  Particulate  materials 

Particulate  material  refers  to  objects  ranging  in  size  from  boulders  in  a 
rock  fill  dam  to  submicron  sized  alloy  powders  in  metallic  components.  In 
principle,  these  materials  are  very  different,  but  the  mechanical  behavior 
due  to  the  granular  nature  is,  in  fact,  very  similar.  The  wide  range  of 
objects  that  is  regarded  as  particulate  material  exhibits  similar  behavior 
when  tested  under  similar  conditions,  despite  differences  in  size,  shape, 
and  composition.  More  important  than  the  size  of  particles,  the 
distribution  of  grain  sizes  greatly  influences  granular  behavior  (Hitcher 
1998).  Granular  media  are  a  subset  of  particulate  materials  that  are  visible 
to  the  naked  eye. 
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The  heterogeneous  characteristics  and  configuration  of  grains  are 
important  mechanical  features.  Macroscopic  behavior  is  determined  by 
the  arrangement  and  interaction  of  discrete  grains  (Oda  and  Iwashita 
1999).  Fabric  describes  the  geometry  of  particle  packing  or  microstructure 
(Bathurst  and  Rothenburg  1988).  The  fabric  of  granular  materials 
provides  insight  on  quantitative  relationships  between  fabric  and 
properties,  mechanics  of  strength  mobilization,  nature  of  peak  and 
residual  strength,  stress-strain  behavior,  and  the  effects  of  different 
sampling  methods  (Kuo  and  Frost  1997). 

Particulate  materials  are  composed  of  mutually  contacting  solid  particles 
that  exhibit  shear  induced  volume  change  and  sensitivity  to  hydrostatic 
stress.  The  behavior  of  particulate  materials  is  controlled  by  the  composi¬ 
tion  and  dynamics  of  the  composing  particles.  Particulate  materials 
behave  in  a  number  of  different  modes  including  deforming  as  a  solid, 
flowing  as  a  fluid,  or  moving  independently  as  individual  particles.  The 
multiple  exhibited  behavior  modes  are  due  to  three-phase  (solid,  liquid, 
gas)  composition,  discontinuous  arrangement,  and  heterogeneous  nature 
(Oda  and  Iwashita  1999).  The  structural  units  of  a  particulate  material  are 
in  mutual  contact,  which  restricts  the  motion  of  each  individual  structural 
unit  (Hitcher  1998).  The  degree  of  motion  restriction  is  controlled  by  the 
magnitude  of  the  applied  stress.  Subjected  to  small  mean  stresses  granular 
materials  flow  like  liquids;  however,  when  subject  to  high  stresses  the 
materials  can  support  large  loadings  (Hitcher  1998).  Extremely  complex 
continuum  models  are  required  to  account  for  all  of  these  different  modes 
and  requiring  sophisticated  model  calibration  techniques. 

2.2  Macro/micro-scale  interaction 

The  micromechanical  behavior  of  particulate  systems  is  commonly 
described  using  macro-scale  measures  of  stress  and  strain.  Historically, 
constitutive  equations  have  been  developed  to  describe  the  behavior  of 
granular  materials.  With  respect  to  soils,  determination  of  the  stress  state 
under  a  particular  loading  configuration  and  proximity  to  the  failure  plane 
in  the  stress  space  are  required  for  the  design  and  evaluation  of  soil 
systems.  The  Mohr-Coulomb  failure  criteria  are  typically  used  to  describe 
the  critical  combination  of  normal  and  shear  stresses  that  define  the 
failure  plane  for  a  given  material.  However,  the  internal  friction  angle  used 
in  the  equation  is  not  a  material  constant,  but  varies  with  void  ratio,  fabric, 
and  stress  state  (Oda  and  Iwashita  1999).  Applications  that  require 
knowledge  of  the  stresses  in  a  soil  system  include  estimates  of 
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compressibility,  bearing  capacity,  stability,  and  lateral  pressure  applied  to 
a  retaining  structure.  The  behavior  of  particulate  materials,  and  in 
particular  soils,  is  complex,  and  the  development  of  a  universally  accepted 
model  has  yet  to  occur  (Jiang  and  Yu  2006).  The  development  of  a 
universal  model  is  obstructed  by  complex  behaviors  observed  in  granular 
materials  including  induced  anisotropy,  phase  transformation,  and  non¬ 
homogeneity  (Sitharam  et.  al  2002).  The  mechanics  associated  with  large- 
scale  soil  deformation  are  poorly  understood  and  traditional  models  are 
not  applicable  (Horner  et.  al  2001). 

The  macroscale  properties  of  an  assembly  of  particles  including  coefficient 
of  contact  friction,  void  ratio,  particle  shape,  and  size  distribution  are  often 
used  to  describe  observed  mechanical  behavior.  The  measurement  of 
mechanical  properties  is  more  difficult  and  costly  than  the  measurement 
of  descriptive  properties  of  particulate  materials  (Feda  1982).  However, 
the  underlying  micromechanics  are  the  true  drivers  of  the  system.  The 
distinguishing  characteristic  of  particulate  material  is  that  deformations 
occur  as  a  result  of  sliding  and  realignment  of  particles  versus  deformation 
of  the  units  themselves.  Therefore,  observed  macro-scale  deformation  is  a 
result  of  sliding  and  rolling  at  particle  contact  points.  When  observed  at 
the  microscale,  the  interaction  of  particles  is  simple  with  explicit  material 
constants  (Oda  and  Iwashita  1999).  Particles  transmit  forces  through  a 
finite  number  of  active  critical  contacts  with  neighboring  particles  (Chris- 
toffersen  et.  al  1981).  The  number  and  strength  of  contacts  determine  the 
strength  and  rigidity  of  a  particulate  material.  Particles  interact  according 
to  contact  force/inter-particle  compliance  relationships  at  the  contacts. 
Despite  the  apparent  simplicity  of  the  system,  the  behavior  is  quite 
complex  at  the  microscopic  level  (Bathurst  and  Rothenburg  1988). 

2.3  Laboratory  methods 

The  shear  strength  of  a  soil  mass  determines  the  resistance  to  failure  and 
sliding  along  the  failure  plane.  Determination  of  strength  parameters  is 
often  accomplished  through  direct  shear  or  triaxial  compression  testing. 
The  direct  shear  test  is  a  simple  and  economical  method  for  determining 
the  strength  parameters  of  soils.  However,  this  method  has  a  number  of 
shortcomings  including  reliability  of  results  and  the  non-uniform 
distribution  of  stresses  (Das  2010).  The  triaxial  shear  test  is  widely  used  in 
both  research  and  conventional  applications.  This  method  is  regarded  as 
highly  reliable,  although  the  equipment  is  costly  and  testing  is  time 
consuming  relative  to  the  direct  shear  test  (Das  2010).  Methods  for 
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conducting  triaxial  testing  include  consolidated-drained  (CD), 
consolidated-undrained  (CU),  and  unconsolidated-undrained  (UU) 
testing.  Consolidated-drained  testing  is  often  used  to  determine  the  long¬ 
term  stability  of  soil  embankments.  The  CU  test  is  used  to  examine  the 
stability  of  consolidated  structures  undergoing  rapid  loading. 
Unconsolidated-undrained  testing  is  used  to  evaluate  the  end-of- 
construction  stability  of  saturated  cohesive  soils. 

2.4  Numerical  methods 

The  primary  advantage  of  numerical  techniques  versus  physical  experi¬ 
ments  is  the  generation  of  abundant  information  on  particle  displace¬ 
ments,  contact  forces,  and  other  physical  quantities.  Numerical  tests  can 
be  used  to  supplement  real  tests  investigating  the  mechanical  behavior  of 
granular  materials  (Oda  and  Iwashita,  1999). 

Continuum  theory  formulates  engineering  problems  as  mathematical 
boundary  value  problems.  A  major  assumption  of  continuum  theory  is  that 
material  behavior  observed  in  small  laboratory  tests  can  be  extrapolated  to 
large  masses  in  the  field.  Hitcher  (1998)  comments  that  the  Mohr-Coulomb 
criterion  is  merely  an  approximation  of  intermediate  stress.  Continuum 
approaches  encounter  issues  due  to  instability  and  limitations  in  the 
constitutive  equations.  Continuum  models  employ  a  number  of  materials 
constants,  some  with  no  clear  physical  meaning  (Oda  and  Iwashita  1999). 
Continuum  mechanics  may  not  be  sufficient  to  express  the  overall 
mechanical  properties  of  granular  materials  (Oda  and  Iwashita  1999). 

Discrete  modeling  can  be  used  to  develop  a  fundamental  understanding  of 
material  behavior  that  can  aid  in  the  development  of  physics  based 
constitutive  models.  Discrete  element  approaches  provide  more  realistic 
models  of  granular  material  behavior  (Oda  and  Iwashita  1999).  The  primary 
limitation  in  applying  discrete  modeling  techniques  to  soil  mechanics  is  the 
number  of  particles  required  to  simulate  problems  of  engineering 
significance.  A  cubic  centimeter  of  clay  contains  roughly  5X1012  particles 
while  the  most  powerful  computing  platforms  running  in  parallel  can 
effectively  handle  around  10?  particles.  Limited  computational  ability  has 
been  overcome  in  previous  studies  by  approximating  the  behavior  of  a 
volume  of  granular  material  using  a  single  spherical  element  (Horner  et.  al 
2001). 


ERDC/GSL  TR-16-29 


8 


Bardet  (1998)  identifies  eight  discrete  numerical  techniques  including 
distinct  element,  modal,  momentum-exchange,  multibody  dynamics, 
structural  mechanics,  mean  field,  and  energy  minimization  methods  in 
addition  to  discontinuous  deformation  analysis.  The  original  DEM,  the 
distinct  element  method,  was  originally  developed  by  Cundall  (1974)  and 
with  minor  modifications,  is  still  the  predominate  DEM  used  today.  The 
distinct  element  method  is  based  upon  the  finite  difference  formulation  of 
the  equations  of  motion.  The  governing  equations  for  the  distinct  element 
method  yield  resultant  forces  and  moments  at  the  center  of  particles 
resulting  from  contact,  body,  inertia,  and  damping  forces.  Assumptions  of 
the  distinct  element  method  include  particles  are  rigid  spherical  bodies, 
contact  points  between  particles  occur  over  an  infinitesimally  small  area, 
particles  can  overlap  slightly  at  contacts,  the  magnitude  of  the  contact 
force  is  determined  by  the  overlap  magnitude,  and  slip  between  particles 
occurs  according  to  Mohr-Coulomb  law  (Oda  and  Iwashita  1999).  DEM 
requires  relatively  few  calibration  parameters  and  gives  realistic  results;  a 
constitutive  equation  requires  many  (>15)  parameters  and  only 
approximates  behavior  under  a  limited  range  of  stress  paths. 

2.5  Particle-size  distribution 

The  behavior  of  a  granular  material  is  strongly  influenced  by  the  grain  size 
distribution  of  the  composing  particles  (Evans  et  al.  2009).  Armstrong  and 
Dunlap  (1966)  comment  that  grain  size  distribution  plays  a  significant  role 
in  particle  interlocking  and  therefore,  will  have  a  definite  effect  on  the 
stress-strain  behavior.  The  primary  effort  in  the  design  of  granular 
material  based  technologies  is  the  mixture  design  process.  Successful 
mixture  design  is  based  on  choosing  proportions  of  coarse  and  fine 
materials  to  optimize  desired  mechanical  properties  (Peters  and  Berney 
2010).  The  influence  of  grain  size  distribution  applies  to  all  particulate 
materials  regardless  of  size.  Additionally,  it  has  been  observed  that  the  size 
ratio  of  particles  dictates  the  packing  density,  not  the  absolute  size 
(McGeary  1961). 

Vallejo  (2001)  conducted  laboratory  and  theoretical  analyses  to  determine 
the  influence  of  grain  size  distribution  on  developed  porosity  in  bimodal 
blends  of  spherical  glass  beads.  The  author  hypothesized  that  the  properties 
and  proportions  of  the  composing  binary  blend  of  particles  determine  the 
macro-scale  observable  behavior.  The  beads  consisted  of  5.0-  and  0.4-mm 
smooth  glass  with  a  specific  gravity  of  2.55.  Laboratory  testing  consisted  of 
uniaxial  compression  and  direct  shear  testing  at  three  different  stress  levels. 
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Findings  of  the  Vallejo  study  included  the  fine  material  in  bimodal  mix¬ 
tures  reduces  point  to  point  contacts  between  coarse  particles,  preventing 
rotation  and  reorientation,  to  stabilize  the  assembly.  Either  the  fine  or 
coarse  material,  depending  on  the  relative  concentrations,  controls  the 
shear  strength  of  bimodal  distributions.  A  blend  with  greater  than 
70  percent  coarse  material  behaves  as  a  uniformly  coarse  material.  Blends 
with  less  than  40  percent  coarse  material  behave  as  a  uniformly  fine 
material.  Blends  between  70  and  40  percent  coarse  material  exhibit 
transitional  behavior  dominated  by  contributions  from  both  the  fine  and 
coarse  material.  The  change  in  behavior  observed  at  the  previous  limits  is 
explained  by  the  porosity  of  the  different  proportions.  Particle  angularity, 
surface  roughness,  and  size  distribution  also  influence  the  porosity  of 
granular  materials.  Eliminating  angularity  and  surface  roughness  through 
the  use  of  smooth  spherical  particles  allows  for  analysis  based  on  size 
distribution  alone.  The  minimum  porosity  condition  is  reached  when  the 
voids  in  the  coarse  matrix  are  completely  filled  with  fine  material.  Binary 
mixtures  can  be  classified  as  coarse  grain  supported,  transitional  coarse- 
fine  grain  supported,  and  fine  grain  supported  structures  depending  on 
the  coarse  material  concentration.  Shear  strength  of  the  binary  mixtures 
mirrored  the  relationship  observed  for  porosity  where  maximum  shear 
strength  was  observed  at  the  same  relative  concentrations  as  minimum 
porosity  (Vallejo  2001). 

Deng  and  Xiao  (2010)  investigated  and  modeled  the  stress-strain  charac¬ 
teristics  of  mixtures  of  expanded  polystyrene  (EPS)  beads  blended  with 
sand  at  varying  ratios  of  EPS  content.  Laboratory  testing  consisted  of 
consolidated-drained  triaxial  compression  testing.  The  EPS-sand  mixtures 
were  blended  and  compacted  into  1.5-in.  by  3.1-in.  specimen.  Specimens 
created  with  varying  EPS  content  were  tested  in  triaxial  compression  at 
four  different  confining  pressures.  The  researchers  observed  a  shift  in  the 
dominant  shearing  mechanism  from  dilative  with  well-defined  peak 
strength  to  fully  contractive  with  no  indicated  peak,  as  the  EPS  bead  con¬ 
tent  was  increased. 

2.6  Percolation  threshold 

The  percolation  threshold  in  particulate  materials  references  the  critical 
proportioning  of  materials  that  allows  for  a  transition  in  mass  material 
behavior.  The  characteristics  of  a  particulate  mass  are  generally  dictated  by 
the  properties  of  the  most  abundant  material.  The  critical  proportioning 
that  causes  the  behavior  of  the  blend  to  be  dominated  by  one,  versus  the 
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other  proportion,  are  determined  by  the  properties  of  the  composing 
materials  (Peters  and  Berney  2010).  The  existence  of  the  percolation 
threshold  derives  from  percolation  theory,  which  describes  behavior 
associated  with  connectivity.  The  percolation  threshold  is  characterized  by 
dramatic  changes  in  system  behavior  resulting  from  small  changes  in 
constituent  composition  fractions  (Peters  and  Berney  2010).  Thornton 
(2000)  theorized  sudden  drops  in  coordination  number  may  also  indicate 
the  percolation  threshold,  resulting  in  a  transition  in  mechanical  behavior. 


Complimentary  to  the  percolation  threshold  is  the  development  and 
evolution  of  force  chains.  With  respect  to  a  binary  granular  material  with 
both  coarse  and  fine  components,  once  the  percolation  threshold  is 
reached,  sufficient  connectivity  exists  within  the  coarse  fraction  for  the 
formation  of  force  chains.  Force  chains  develop  in  the  large  fraction  as  the 
proportion  becomes  significant  enough  for  forces  to  be  transmitted 
exclusively  through  the  large  fraction  of  the  medium.  Figure  1  illustrates  a 
bimodal  blend  with  behavior  dominated  by  fine  proportion  (A),  increasing 
the  connectivity  of  the  coarse  proportion  (B),  and  sufficient  coarse 
proportion  for  the  formation  of  force  chains  and  shift  in  dominant 
proportion  (C).  When  the  percolation  threshold  of  materials  such  as  sand 
and  silty-clay  blends  as  studied  by  Peters  and  Berney  (2010)  is  reached, 
erratic  and  even  unstable  behavior  was  observed. 


Figure  1.  Influence  of  particle  proportioning  on 
stress  transfer  mechanisms. 
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3  Discrete  Element  Method 

Inter-particle  forces  greatly  influence  granular  material  behavior,  but  it  is 
virtually  impossible  to  measure  inter-particle  forces  in  real  granular 
materials  (Peters  and  Berney  2010).  Due  to  the  difficulty  in  directly 
measuring,  the  forces  are  typically  estimated  from  the  boundary  forces 
measured  in  physical  testing.  The  discrete  element  method  (DEM)  is  a 
numerical  method  for  simulating  the  behavior  of  particle-based  systems 
by  integrating  the  equations  of  motion  for  each  body  in  an  assemblage  of 
particles  (Cundall  and  Strack  1979).  The  method  treats  each  particle 
within  the  assemblage  as  a  rigid,  discrete  body.  Binary  contact  laws  define 
the  forces  and  motions  that  occur  in  the  system.  Contact  forces  and 
displacements  are  calculated  for  equilibrium  conditions  where 
incremental  changes  in  contact  force  are  determined  by  incremental 
displacements  at  the  particle  contact.  Newton’s  second  law  of  motion  is 
integrated  to  obtain  the  net  forces  and  moments  acting  on  each  particle. 
This  approach  deviates  from  other  numerical  methods  such  as  finite 
elements  in  that  the  observed  macro-scale  behavior  is  a  result  of  the 
definition  of  particle  interactions  versus  a  definition  of  the  assemblage. 

DEM  refers  to  a  number  of  numerical  computational  methods  for 
modeling  granular  materials.  Bardet  (1998)  describes  eight  and  O’Sullivan 
(2002)  documented  15  different  discrete  numerical  methods.  Advantages 
of  DEM  include  the  ability  to  model  a  wide  variety  of  complex  granular 
flow  and  rock  mechanics  applications  while  generating  detailed  data  on 
the  micro-mechanics  occurring  at  the  particle  scale  not  obtainable  using 
physical  experiments.  The  primary  disadvantage  of  DEM  is  that  the 
method  is  computationally  intensive  and  thus,  scaled  models  of  actual 
systems  incorporating  thousands  of  particles  versus  billions  must  be  used 
in  order  to  obtain  results  in  a  reasonable  time  period.  Cundall  (1974) 
developed  the  numerical  method  known  as  the  distinct  element  method  to 
describe  the  mechanical  behavior  of  granular  materials.  The  distinct 
element  method  is  an  explicit  scheme  in  which  the  behavior  of  a  system  of 
particles  is  determined  by  interactions  at  shared  contacts. 

The  distinct  element  method  is  designed  for  applications  where  large 
displacements  are  of  interest.  The  primary  assumption  of  the  distinct 
element  method  is  that  macro-scale  deformations  are  a  result  of  the  rolling 
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and  sliding  of  spherical  particles  and  not  deformation  of  the  particles 
themselves.  The  assumption  of  non-deformable  particles  is  valid  in  low- 
stress  situations  where  rotation  and  sliding  are  the  primary  mechanisms  of 
movement.  However,  in  conditions  where  the  magnitude  of  particle 
loading  is  very  high  inducing  particle  deformation  of  a  similar  magnitude 
to  the  deformation  in  the  system,  the  distinct  element  method  models  will 
not  be  valid  and  another  approach,  such  as  the  finite  element  method, 
should  be  investigated  (Cundall  1974).  Cundall  and  Strack  (1979)  observed 
that  precise  modeling  of  particle  deformation  was  not  required  to 
accurately  simulate  mechanical  behavior.  An  additional  assumption  of  the 
method  is  that  particles  interact  at  the  point  of  contact  and  forces  are 
adjusted  until  a  state  of  equilibrium  is  reached.  Cundall  (1974)  comments 
that  the  assumption  that  particle  interactions  can  be  adequately  modeled 
by  two  localized  forces  seems  more  valid  than  the  conventional  assump¬ 
tion  of  a  continuous  uniform  stress  because  real  particles  typically  only 
touch  at  a  few  asperities  regardless  of  surface  smoothness. 

3.1  DEM  in  geomechanics  research 

The  DEM  has  been  successfully  used  by  a  number  of  researchers  to  provide 
insight  on  the  micromechanics  of  granular  systems  while  accurately 
simulating  macro-scale  behavior.  A  summary  of  select  works  documenting 
the  utilization  of  DEM  in  geomechanics  research  is  presented  in  Table  1. 

The  DEM  was  originally  developed  as  a  two-dimensional  (2-D)  numerical 
method  for  predicting  deformations  in  rock  structures  (Cundall  1974). 

Other  researchers  including  Cundall  and  Strack  (1979),  Jiang  and  Yu 
(2006),  and  Evans  et  al.  (2009)  have  applied  2-D  DEM  for  deeper  under¬ 
standing  of  the  underlying  micromechanics  occurring  in  granular  materials. 

Cundall  and  Strack  (1979)  validated  the  performance  of  discrete  numerical 
models  by  reproducing  historical  results  from  a  photoelastic  analysis  of  an 
assemblage  of  discs.  The  photoelastic  analysis  consisted  of  loading  two 
horizontal  and  two  vertical  beams  that  formed  the  boundaries  of  an 
assemblage  of  discs.  The  bounding  beams  were  loaded  and  allowed  to  freely 
rotate.  The  motion  of  the  discs  was  recorded  to  construct  force  vector  plots. 
The  physical  experiment  was  simulated  using  the  distinct  element  method 
to  obtain  a  qualitative  assessment  of  the  effectiveness  of  discrete  numerical 
models  for  modeling  the  behavior  of  particulate  systems.  The  authors 
observed  good  correspondence  between  the  physical  test  and  numerical 
simulation  and  concluded  the  distinct  element  method  is  effective  for 
modeling  the  fundamental  behavior  of  granular  assemblies  (Cundall  and 
Strack  1979). 
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Table  1.  Summary  of  select  geomechanics  works  employing  DEM. 


Author(s) 

Title 

Keywords 

P.A. 

Cundall 

A  Discrete  Numerical  Model  for 

2-D  DEM,  discs,  axisemmetric 

O.D.L. 

Strack 

Granular  Assemblies 

compression,  photoelastic  analysis 

M. 

Jiang 

Application  of  Discrete  Element 

2-D  DEM,  discs,  grain  size  distribution, 

H.S. 

Yu 

noncoaxality,  effective  stress,  bonding. 

penetration, 

T.  Matthew 

Evans 

Hamed 

Mojarrad 

Grain  Size  Distribution  Effects  in  2D 

2-D  DEM,  discs,  biaxial  compression, 

Charles 

Cunningham 

Discrete  Numerical  Experiments 

grain  size  distribution 

Akhtar  A. 

Tayebali 

C. 

Thornton 

3-D  DEM,  spherical  particles, 

Numerical  Simulations  of  Deviatoric 

axisemmetric  compression,  packing 

Shear  Deformation  of  Granular  Media 

density,  contact  force,  coordination 

number,  anisotropy 

S.  Joseph 

Anthony 

3-D  DEM,  spherical  particles,  cubical 
specimen,  axisemmetric  compression, 
packing  density,  interface  energy 

Evolution  of  Force  Distribution  in 

Three-Dimensional  Granular  Media 

David  A. 

Homer 

3-D  DEM,  spherical  particles,  finite 

John  F. 

Peters 

Large  Scale  Discrete  Element  Modeling 

elements,  high  perfonnance  computing. 

AlexR. 

Carrillo 

of  Vehicle-Soil  Interaction 

parallelization,  particle  interaction. 

mass  deformation 

L. 

Cui 

Exploring  the  Macro-  and  Micro-Scale 

3-D  DEM,  spherical  particles,  direct 
shear,  steel  spheres 

C. 

O'Sullivan 

Response  of  an  Idealized  Granular 

Material  in  the  Direct  Shear  Apparatus 

Thallak  G. 

Sitharam 

Micro  mechanical  Modeling  of 

3-D  DEM,  spherical  particles,  triaxial 

S.V. 

Dinesh 

Monotonic  Drained  and  Undrained 

compression,  stress  path,  packing 

N. 

Shimizu 

Shear  Behavior  of  Granular  Media 

density,  coordination  number, 

using  Three-Dimensional  DEM 

anisotropy 

Catherine 

O'Sullivan 

Examination  of  the  Response  of 

Jonathan  D. 

Bray 

Regularly  Packed  Specimens  of 

3-D  DEM,  spherical  particles,  triaxial 

Michael 

Riemer 

Spherical  Particles  Using  Physical 

compression,  uniform  and  non-uniform 
size  distribution,  steel  spheres 

Simulations 

Catherine 

O'Sullivan 

2-D/3-D  DEM,  circular  rods,  spherical 

Jonathan  D. 

Bray 

Experimental  Validation  of  Particle- 

particles,  non-uniform  size  distribution, 

Liang 

Cui 

Based  Discrete  Element  Methods 

biaxial  compression,  triaxial 

compression,  direct  shear 

L. 

Cui 

An  Analysis  of  the  Triaxial  Apparatus 

3-D  DEM,  spherical  particles,  triaxial 

C. 

O'Sullivan 

using  a  Mixed  Boundary  Three- 

compression,  steel  spheres,  contact 

s. 

O'Neill 

Dimensional  Discrete  Element  Model 

force,  fabric 

Tang-Tat 

Ng 

Discrete  Element  Simulations  of  the 

Critical  State  of  a  Gtanular  Material 

3-D  DEM,  ellipsoidal  particles,  cubical 
specimen,  triaxial  compression 

N. 

Belheine 

J.P. 

Plassiard 

Numerical  Simulation  of  Drained 

3-D  DEM,  spherical  particles,  triaxial 
compression,  Labenne  sand 

F.V. 

Donze 

Triaxial  Test  using  3D  Discrete  Element 

F. 

Darve 

Modeling 

A.. 

Seridi 
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Jiang  and  Yu  (2006)  present  applications  of  DEM  to  geomechanics.  Topics 
covered  are  based  on  the  authors’  previous  research  and  include  nonco¬ 
axiality  of  granular  materials,  effective  stress  in  unsaturated  soils,  bonding 
effect  in  natural  soils,  and  penetration  mechanisms  in  granular  ground. 
The  DEM  model  used  in  the  study  was  a  modified  version  of  the  code  pre¬ 
sented  by  Cundall  and  Strack  (1979).  The  authors  observed  good  agree¬ 
ment  between  DEM  simulations  and  laboratory  simple  shear  testing.  The 
authors  used  the  DEM  as  a  virtual  laboratory  to  gather  data  not  possible  in 
physical  testing.  The  DEM  was  observed  to  be  an  effective  tool  for  examin¬ 
ing  the  noncoaxiality  of  granular  materials,  effective  stress  in  unsaturated 
granular  materials,  bonding  effect  in  natural  soils,  and  complex  boundary 
value  problems  in  geotechnical  engineering. 

Evans  et  al.  (2009)  investigated  the  influence  of  grain  size  distribution  on 
the  stress-strain-strength  behavior  of  2-D  systems  of  particles.  A  DEM 
numerical  model  was  used  to  simulate  biaxial  compression  testing  on  four 
distinct  grain  size  distributions  including  coarse,  medium,  and  fine  well- 
graded  and  uniformly  graded  distributions.  Elliptical  particles  with  an 
aspect  ratio  of  1.5  to  1.0  were  used  in  the  model  and  were  composed  by 
bonding  pairs  of  identical  spheres.  The  authors  concluded  that  DEM  is 
capable  of  modeling  the  response  of  well-graded  assemblies  in  addition  to 
the  more  simple  gradations  typically  used. 

The  accurate  simulation  of  physical  systems  requires  accurate  representa¬ 
tion  of  the  modeled  media  in  three-dimensional  (3-D).  A  number  of 
researchers  have  extended  successful  2-D  simulations  of  granular  material 
behavior  into  3-D.  Thornton  (2000),  Antony  (2000),  Horner  et  al.  (2001), 
and  Cui  and  O’Sullivan  (2006)  have  all  used  3-D  DEM  to  investigate  the 
behavior  of  granular  media. 

Thornton  (2000)  investigated  the  use  of  DEM  models  to  simulate  3-D 
polydisperse  systems  of  elastic  spheres  in  axisymmetric  compression, 
radial  deviatoric  loading,  constant  deviatoric  strain,  and  multi-axial  plane 
strain  testing.  The  specimens  consisted  of  five  discrete  sized  particle  dis¬ 
tributions.  The  author  observed  the  stress-strain-dilation  response  of 
simulations  was  similar  to  laboratory  observations  of  axisymmetric  com¬ 
pression.  Thornton  concluded,  “Realistic  3-D  stress-strain  dilation  behav¬ 
ior  can  be  generated  by  discrete  element  simulations.” 
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Antony  (2000)  investigated  the  evolution  of  contact  normal  force  distri¬ 
bution  in  monodisperse  particulate  systems  using  DEM  simulations.  The 
DEM  model  included  a  submerged  particle  system  containing  5,000  uni¬ 
form  spheres  in  a  cubical  periodic  cell  undergoing  axisymmetric  compres¬ 
sion.  Variables  included  noncohesive  particles,  low  interface  energy, 
packing  density,  and  stiffness.  The  transmission  of  forces  during  slow 
shearing  was  analyzed.  The  author  found  DEM  to  be  an  effective  tool  for 
investigating  the  micromechanics  of  particulate  systems. 

Horner  et.  al  (2001)  present  the  promise  and  associated  concerns  when 
parallelized  DEM  modeling  is  applied  to  discontinuous,  large  deformation 
soil-vehicle  interaction  problems.  The  authors  modeled  the  deformation  in 
a  soil  mass  resulting  from  the  interaction  with  a  tined  plow  using  a  3-D 
DEM  consisting  of  10,000,000  particles  optimized  for  use  on  a 
parallelized  high  performance  computing  platform.  The  authors  comment 
on  the  value  of  modeling  vehicle  design  alternatives  in  a  virtual 
environment  and  on  the  difficulty  of  using  continuum  based  methods.  The 
particulate  nature  of  the  DEM  allows  for  modeling  discontinuous  large 
deformations  encountered  in  such  problems.  The  authors  applied  a  joint 
DEM  finite  element  method  (FEM)  model  to  accurately  simulate  the 
phenomenon  observed  in  field-testing.  The  authors  concluded  that  DEM 
captures  the  particulate  phenomenology  of  soils  and  avoids  the 
programming  complexity  of  traditional  constitutive  methods.  The 
researchers  also  observed  a  remarkable  level  of  detail  in  DEM  simulations 
of  large-scale  deformation  also  observed  in  physical  testing  (Figure  2). 
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Figure  2.  DEM  of  tined  plow  moving  through  particle  system  with  colorization  based 
on  stress  intensity  (Horner  et.al.  2001). 


Cui  and  O’Sullivan  (2006)  examined  the  similarities  and  differences  in 
micro-  and  macro-performance  of  idealized  granular  assemblies  composed 
of  uniform  steel  spheres.  Laboratory  testing  included  direct  shear  tests  on 
steel  sphere  assemblies  with  varying  normal  load  at  a  fixed  velocity.  DEM 
simulations  of  the  direct  shear  tests  were  attempted  using  material  prop¬ 
erties  determined  in  the  laboratory.  The  macro-performance  of  the  phys¬ 
ical  and  numerical  testing  were  compared,  and  micro-mechanical 
performance  measures  including  distribution  of  force,  displacements, 
rotations,  specimen  fabric,  and  variation  in  volumetric  shear  strain  deter¬ 
mined  using  the  DEM  model.  The  DEM  model  was  found  to  effectively 
simulate  the  direct  shear  test.  The  model  did  under-predict  the  peak 
shearing  resistance  angle;  however,  the  authors  hypothesized  that  the  dif¬ 
ferences  may  be  caused  by  non-constant  inter-particle  friction  in  physical 
models  (increasing  with  strain).  The  DEM  confirmed  theories  about 
inhomogeneities  in  stress,  strain,  and  rotation  in  the  direct  shear  test. 

Triaxial  compression  testing  is  the  primary  method  for  determining 
strength  and  deformation  behavior  in  geomechanics.  Most  commonly 
accepted  constitutive  theories  for  soil  behavior  have  been  developed  using 
extensive  triaxial  testing.  A  number  of  researchers  including  Sitharamet.  al 
(2002),  O’Sullivan  et.  al  (2004),  O’Sullivan  et.  al  (2006),  Cui  et.  al  (2007), 
Ng  (2009),  and  Belheine  et  al.  (2009)  have  used  3-D  DEM  simulations  of 
triaxial  testing  to  develop  an  understanding  of  the  underlying  microme¬ 
chanics  of  soils  subjected  to  external  loading. 
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Sitharam  et.  al  (2002)  investigated  the  evolution  of  induced  anisotropy 
during  drained  and  undrained  triaxial  tests  on  loose  and  dense  assemblies 
of  spheres  using  a  DEM  model.  Simulations  were  conducted  on  assemblies 
of  1,000  spherical  particles  with  lognormal  particle-size  distributions 
representing  Unified  Soil  Classification  System  (USCS)  GP  and  SP  soils. 
Simulation  of  the  triaxial  testing  of  cubical  specimens  was  accomplished 
with  a  focus  on  the  evolution  of  microscopic  features  of  the  granular 
assembly.  The  simulation  of  triaxial  undrained  tests  was  accomplished 
with  DEM  by  maintaining  a  constant  volume  of  the  specimen  during 
shear.  The  authors  comment  on  the  difficulty  in  modeling  the  complex 
behavior  of  soils  including  induced  anisotropy,  phase-transformation,  low- 
confinement  behavior,  and  non-homogeneity  using  traditional  continuum 
mechanics  approaches.  In  granular  materials,  such  as  soils  that  are 
dominated  by  particulate  behavior,  it  is  advantageous  to  treat  the  medium 
as  an  assemblage  of  particles  rather  than  a  continuum.  The  authors 
conclude  that  DEM  models  are  valuable  in  geomechanics  because  they 
provide  complete  quantitative  data  on  all  microscopic  features  of  an 
assembly  of  particles. 

O’Sullivan  et.  al  (2004)  investigated  the  response  of  spherical,  idealized 
granular  assemblies  in  both  physical  and  numerical  experiments.  Physical 
testing  included  plain  strain  and  triaxial  compression  tests.  The  numerical 
experiments  employed  a  DEM  model  to  simulate  the  physical  testing  by 
matching  the  properties  of  the  spheres  observed  in  physical  testing  to 
include  friction  angle,  fabric  (void  ratio,  coordination  number),  and 
packing  configuration.  The  DEM  simulations  were  able  to  accurately 
capture  the  peak  response  of  the  tested  materials.  However,  the  post-peak 
behavior  was  difficult  to  model  and  was  determined  to  be  a  function  of 
discrepancies  in  specimen  fabric.  The  researchers  concluded  that  the  DEM 
was  effective  for  monitoring  both  changes  in  coordination  number  and 
strength,  with  the  strength  of  the  specimen  decreasing  with  decreasing 
coordination  number. 

O’Sullivan  et.  al  (2006)  set  out  to  quantitatively  demonstrate  that  DEM 
models  could  accurately  simulate  the  response  of  granular  materials.  A 
series  of  three  DEM  validation  tests  were  undertaken  in  the  investigation. 
The  first  study  examined  the  2-D  response  of  uniform  and  almost  uniform 
hexagonally  packed  rods.  The  biaxial  compression  of  2-D  hexagonally 
packed  rods  was  examined  in  the  laboratory  and  modeled  using  DEM.  The 
second  study  was  an  extension  of  the  first  looking  at  the  3-D  response  of 
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uniform  spheres  in  face-centered-cubic  (FCC)  and  rhombic  packing 
arrangements  in  triaxial  and  plane  strain  compression.  The  culminating 
study  examined  the  response  of  randomly  packed  spheres  in  the  direct 
shear  test.  The  DEM  models  were  found  to  be  able  to  accurately  replicate 
the  response  of  granular  materials  in  standard  laboratory  tests  while 
providing  insight  into  the  particle-scale  interactions  driving  the  observed 
macro-scale  response.  The  DEM  was  found  to  satisfactorily  simulate 
laboratory  response  up  to  the  peak  strength.  However,  post-peak  response 
could  not  be  captured  without  modification  of  model  parameters.  The 
authors  concluded  DEM  modeling  is  a  valuable  tool  for  engineers 
examining  the  micro-scale  mechanics  driving  the  response  of  granular 
materials. 

Cui  et.  al  (2007)  studied  the  micro-scale  response  of  systems  of  steel 
spheres  in  triaxial  testing  using  both  physical  testing  and  DEM  modeling. 
Mixed  boundary  3-D  DEM  code  was  used  in  the  study  to  accurately 
simulate  the  latex  membrane  and  imposed  confining  pressure  used  in  the 
physical  triaxial  test.  Testing  was  performed  on  assemblies  of  uniform  and 
tri-modal  discrete  sphere  sizes  in  face-centered-cubic  (FCC)  packing 
arrangement.  Coupling  DEM  simulations  and  triaxial  tests  on  ideal 
spheres  allowed  for  accurate  replication  of  the  material,  facilitating 
quantitative  validation  of  the  model.  The  authors  commented  that  the 
quality  of  the  validation  would  have  been  compromised  by  the  use  of  more 
complex  materials  such  as  real  soils.  It  was  recommended  to  validate  DEM 
models  using  simple  materials  prior  to  incorporating  more  complex  “real” 
materials  to  ensure  accuracy  and  appropriateness  of  the  model. 

Ng  (2009)  uses  a  series  of  DEM  simulations  to  examine  the  critical  state  of 
three  cubical  assemblies  of  ellipsoids.  The  critical  state  was  reached  by 
simulating  axial  compression,  lateral  extension,  and  constant  mean-stress 
compression  tests  at  varying  confining  pressures.  The  critical  state  line  in 
void  ratio,  mean  stress,  and  deviator  stress  space  was  explored.  The 
numerical  experiments  challenged  conventional  theory  on  the  critical  state 
line.  Curvature  was  observed  in  the  critical  state  line  in  simulations  of 
rigid  materials,  suggesting  that  the  curvature  in  physical  experiments  is 
not  solely  due  to  particle  breakage.  The  author  found  DEM  to  be  an 
effective  tool  for  examining  the  micromechanics  of  soil  structures. 

Belheine  et  al.  (2009)  investigated  the  effectiveness  of  a  DEM  model  in 
reproducing  macro-scale  behavior  of  a  granular  material.  A  DEM  model 
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was  used  to  simulate  drained  triaxial  testing  of  Labenne  sand.  Ten 
thousand  spherical  particles  were  used  in  the  3-D  DEM  model.  The 
authors  found  the  discontinuous  nature  of  granular  materials  to  induce 
unique  behaviors  such  as  anisotropy,  micro -fracture,  and  local  instability. 
The  DEM  model  used  in  the  study  reproduced  the  non-linear  stress-strain 
behavior  of  Labenne  sand.  The  DEM  was  also  able  to  accurately  produce 
the  Mohr-Coulomb  criterion  for  the  modeled  material.  The  researchers 
found  the  numerical  and  experimental  results  to  be  in  good  agreement. 

3.2  Governing  equations 

The  DEM  is  an  iterative,  cyclic  calculation  numerical  model.  The  inter¬ 
particle  contact  forces  are  calculated  for  each  particle  using  the  relative 
velocities  and  the  employed  force-displacement  laws.  The  behavior  at  the 
contact  is  critical  to  effective  modeling  because  the  interaction  of  particles 
is  characteristic  to  the  material  of  interest  (Hitcher  1998).  The  force  and 
moment  sums  are  used  to  compute  the  linear  and  rotational  accelerations 
of  each  particle,  which  are  numerically  integrated  to  yield  velocity,  and 
integrated  again  to  give  displacement  according  to  Newton’s  second  law  of 
motion.  Contact  checks  are  then  performed  prior  to  initiation  of  the  next 
calculation  cycle.  Although  the  DEM  system  is  dynamic,  it  can  be 
approximated  as  static,  if  the  loading  rates  at  the  boundaries  are  small 
enough  that  the  inertial  forces  are  small  compared  to  the  contact  forces. 

Forces  in  the  system  arise  from  deformations  occurring  between  particles, 
assumed  to  occur  at  the  shared  contacts.  The  contact  interaction  of  the 
particles  can  be  modeled  using  a  set  of  springs,  dashpots,  and  a  frictional 
slider  as  shown  in  Figure  3.  Elastic  force/displacement  laws  are  used  in 
the  DEM  and  are  incremental,  where  a  change  in  displacement  produces  a 
change  in  force,  which  is  added  to  the  existing  force  (Cundall  1974).  For 
every  time-step,  the  incremental  normal  and  shear  displacements  for  a 
given  contact  are  evaluated  from  the  incremental  rotational  displacement. 
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Figure  3.  Linear-elastic  particle  contact  model. 


The  normal  and  shear  forces  are  calculated  using  the  existing  and  incre¬ 
mental  forces  according  to  the  force/displacement  equations.  The  incre¬ 
mental  contact  forces  are  resolved  in  x  and  y  forces  and  moments  and 
added  to  the  existing  x  and  y  forces  and  moments.  Linear  elasticity  is  suffi¬ 
cient  to  describe  the  behavior  of  grains  except  were  high  stresses  result  in 
plastic  deformation  or  rupture  of  the  grains  (Hitcher  1998).  Elastic 
spheres  interact  in  a  non-linear  manner,  but  a  linear  spring  model  is 
useful  for  verification  of  theoretical  concepts  (Bathurst  and  Rothenburg 
1988). 

To  prevent  excessive  energy  buildup,  a  spring-dashpot  system  is  employed 
to  model  the  conversion  of  kinetic  to  thermal  energy  taking  place  at  the 
contact  points  in  a  physical  system.  The  inclusion  of  damping  allows  the 
assembly  of  particles  to  reach  equilibrium.  Damping  is  achieved  mathe¬ 
matically  by  connecting  viscous  dashpots  in  parallel  with  the  normal  and 
shear  stiffness  at  each  contact  (Cundall  1974).  The  particle  material  is  elas¬ 
tic,  but  through  the  use  of  energy  dissipating  mechanisms  at  the  contact, 
the  overall  behavior  is  inelastic  (Sitharam  et.  al  2002).  No  explicit 
rotational  damping  is  employed  in  the  DEM.  Cundall  and  Strack  (1979) 
commented  that  damping  does  not  influence  the  force  equilibrium  value 
but  reduces  the  number  of  calculation  cycles  required  to  reach  equi¬ 
librium.  Damping  is  incorporated  to  allow  a  state  of  equilibrium  to  be 
reached  in  the  model.  However,  the  departure  from  equilibrium  can  also 
be  minimized  by  reducing  the  loading  rate  (Cundall  and  Strack  1979). 
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3.2.1  Calculation  cycle 

The  calculation  cycle  begins  by  determining  the  relative  velocity  and  rota¬ 
tional  rate  of  each  particle  in  the  assembly  according  to  Equations  l  and  2, 
respectively. 


m 


<9v; 

dt 


mgnf+i^K 

2=1 


(1) 


where: 


vi  =  velocity 
t  =  time  increment 
m  =  mass 

gnf  =  acceleration  due  to  gravity 

Nc  =  number  of  contacts 
f?  =  force  at  the  contact 


*.p 

UL  c= 1  c=l 


where: 


(2) 


Wi  =  rotational  rate 
Im  =  moment  of  inertia 
p  =  unit  weight 

rf  =  radius  from  the  particle  center  to  the  point  of  contact 
mf  =  moment  applied  at  the  contact 


A  determination  of  particle  contact  is  then  conducted  by  evaluating  each 
particle  and  its  nearest  neighbors.  Two  particles  are  in  contact  if  the  dis¬ 
tance  between  the  particle  centers  is  less  than  the  sum  of  the  combined 
radius  as  shown  in  Equation  3. 


d<RA+RB 


(3) 


where: 

d  =  distance  between  the  centers  of  particles  A  and  B 
Ra  =  radius  of  particle  A 
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Rb  =  radius  of  particle  B 


For  two  particles  in  contact,  the  overlap  or  displacement  at  the  contact  is 
determined  according  to  Equation  4. 


=vf-vf 


em 


r?wAk 


(4) 


where: 

Ac{  =  relative  motion  at  the  contact 
vf  =  velocity  of  particle  A 
vf  =  velocity  of  particle  B 

rf  =  radius  from  the  center  of  particle  A  to  the  point  of  contact 
wf  =  rotational  rate  of  particle  A 

rf  =  radius  from  the  center  of  particle  B  to  the  point  of  contact 
wf  =  rotational  rate  of  particle  B 

The  rotational  difference  between  two  particles  in  contact  is  determined  as 
shown  in  Equation  5. 


Aw-  =At(wf  -wf'j 


(5) 


where: 

Aw-  =  difference  in  rotation  of  particles  A  and  B 

At  =  change  in  time  increment 
wf  =  rotational  rate  of  particle  A 
wf  =  rotational  rate  of  particle  B 


The  force-displacement  laws  are  applied  to  determine  the  forces  and 
moments  acting  on  each  particle  as  a  result  of  interaction  with  another 
particle.  The  normal  force,  shear  force,  and  moment  at  each  contact  are 
determined  in  accordance  with  Equations  6  through  8,  respectively. 


/"  = 


Kn  A" 

ErKn[ A°-A"),A"  <A° 


(6) 
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where: 


fn  =  normal  component  of  the  contact  force 
Kn  =  normal  stiffness  constant 

A"  =  normal  component  of  the  contact  displacement 
Er  =  an  energy  dissipation  factor 

A0  =  greatest  magnitude  of  penetration  in  the  history  of  the  normal 
contact  displacement 


ft 


KSA ) 

fn  tan  (pnf,|y;s  |  >  fn  tan  cp 


(7) 


where: 

f.s  =  shear  component  of  the  contact  force 

Ks  =  shear  stiffness  constant 
A*  =  shear  component  of  the  contact  displacement 
cp  =  a  force  friction  parameter 
n?  =  shear  force  unit  vector 


= 


KmAwct 

fn  tan  cp mn;m ,  I m'  >/"  tan  cpn 


(8) 


where: 

m.  =  moment  at  the  contact 

Km  =  moment  stiffness  constant 
n;m  =  moment  unit  vector 


3.2.2  Velocity  verlet  algorithm 

The  velocity  verlet  algorithm  is  a  numerical  method  for  integrating  New¬ 
ton’s  second  law  of  motion.  The  algorithm  is  a  “leapfrog”  approach  that 
uses  the  current  position,  velocity,  and  acceleration  to  calculate  the  future 
state.  The  linear  and  rotational  velocities  are  integrated  to  yield  the  future 
linear  and  rotational  accelerations  and  applied  as  shown  in  Equations  g 
and  to  for  the  future  position  and  velocity  and  Equations  n  and  12  for  the 
angle  of  rotation  and  angular  velocity. 
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rn+l=rn+^n+^t2an 

(9) 

where: 

rn 

=  the  position  of  the  particle  at  the  nth  time-step 

vn 

=  the  velocity 

an 

=  the  acceleration 

Vn+i=vn+^At(an+an+1) 

(10) 

(11) 

where: 

9n  =  the  angle  of  rotation  of  the  particle  at  the  nth  time-step 
wn  =  the  rotational  rate 
an  =  the  angular  acceleration 

^„+i  =  “h+!Af(a„+a„+i)  (12) 

Once  the  position,  velocity,  and  acceleration  of  each  particle  in  the  assem¬ 
bly  for  the  next  time  increment  have  been  determined,  the  calculation 
cycle  is  restarted. 

3.2.3  Contact  detection 

Simplified  DEM  models  achieve  contact  detection  by  comparing  the  dis¬ 
tance  between  particles  with  every  other  particle  in  the  assemblage.  This 
procedure  requires  significant  computational  effort  as  the  number  of  cal¬ 
culations  increases  N2  as  the  total  number  of  particles  (N)  increases.  Most 
practical  DEM  models  employ  a  nearest  neighbor  search  procedure  to 
define  contact  between  particles.  The  nearest  neighbor  approach  assigns 
each  particle  to  a  cell  space  based  upon  its  Cartesian  coordinates.  Thus, 
when  checking  for  contacts  only  the  adjoining  cells  need  be  evaluated.  This 
approach  is  very  simple  for  uniformly  sized  particles;  however,  when  a 
distribution  of  particle  sizes  is  used,  a  bounding  box  cell  search  scheme  is 
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employed.  In  this  method  the  cell  size  is  determined  by  the  smallest 
particle  size,  and  the  bounding  box  is  the  collection  of  all  the  cells  in  which 
any  portion  of  the  particle  crosses  (Peters  et  al.  2010*).  The  contact 
detection  is  conducted  only  on  the  cells  within  the  bounding  box  for  each 
particle  to  ensure  maximum  computational  performance. 

3.2.4  Parameter  selection 

The  material  properties  of  a  continuous  solid  cannot  be  used  directly  when 
selecting  parameters  for  the  DEM  model.  Armstrong  and  Dunlap  (1966) 
propose  the  analogy  that  if  a  steel  cylinder  were  machined  into  spheres, 
the  mass  of  spheres  would  not  exhibit  the  same  modulus  of  elasticity  or 
Poison’s  ratio  as  the  original  cylinder.  The  selection  of  material  properties 
is  accomplished  such  that  the  observed  performance  resembles  physical 
observations.  However,  parameters  should  also  be  selected  to  minimize 
particle  overlapping  in  relation  to  particle  size,  ensure  stability  of  the  sys¬ 
tem,  and  the  values  calculated  in  the  program  are  in  a  range  that  can  be 
handled  with  accuracy  (Cundall  and  Strack  1979).  Cundall  and  Strack 
(1979)  observed  that  while  other  parameters  can  be  adjusted  in  an  effort  to 
optimize  the  program,  the  inter-particle  friction  angle  and  the  ratio  of 
shear  and  normal  stiffness  should  be  selected  based  upon  the  physical 
properties  of  the  material  being  modeled.  O’Sullivan  et.  al  (2004) 
emphasize  the  importance  of  using  the  physical  characteristics  of  the 
material  being  modeled  as  parameters  in  DEM  simulations. 

3.2.5  Critical  time-step 

Failure  in  particulate  systems  is  often  prefaced  by  a  number  of  particle 
realignments  as  the  loads  are  redistributed.  Cundall  (1974)  commented 
that  the  sequence  and  magnitude  of  the  motions  are  important  to  the  con¬ 
ditions  precipitating  final  failure  due  to  the  high  degree  of  non-linearity  in 
a  particulate  system.  Thus,  the  DEM  incorporates  an  explicit  numerical 
method,  where  the  entire  system  of  equations  is  solved  with  small 
increments  of  time,  or  time-steps.  Use  of  the  explicit  method  allows  for 
tracking  the  path  to  failure,  versus  an  implicit  method,  which  is  more 
efficient  from  a  time  perspective  but  may  skip  over  some  of  the  important 
effects. 


*  Peters,  J.  F.,  R.  Kala,  R.  S.  Maier,  L.  Walizer,  J.  Wibowo,  and  R.  E.  Wahl.  2010.  User  guide  for  DEM 
development. 
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The  discrete  element  model  relies  upon  explicit  integration,  which  is  only 
conditionally  stable.  The  DEM  is  based  upon  the  idea  that  a  time  incre¬ 
ment  so  small  can  be  selected  in  the  explicit  scheme  so  that  during  a  single 
increment  interactions  caused  by  a  given  particle  cannot  propagate  beyond 
the  nearest  neighbors  (Cundall  and  Strack  1979).  The  time-step  required 
for  stability,  or  critical  time-step,  is  a  function  of  the  mass  and  stiffness  of 
the  particles,  with  the  most  critical  relationship  as  defined  by  the  natural 
oscillating  frequency  (Cundall  and  Strack  1979).  The  DEM  program  used 
in  this  study  selects  a  time-step  0.10  times  the  critical  time-step,  calculated 
according  to  Equation  13,  and  based  on  the  smallest  and  stiffest  particles 
in  the  system.  The  critical  time-step  limitation  often  requires  compromises 
during  the  calibration  of  the  model  (Peters  et  al.  2010).  To  achieve 
acceptable  run  times  the  stiffness  must  be  assigned  low  or  the  mass 
assigned  high.  Acceptable  results  have  been  achieved  when  applying  either 
approach.  Scaling  particle  mass  does  not  affect  stresses  and  strains  as 
opposed  to  scaling  particle  stiffness.  However,  the  velocities  and 
accelerations  are  reduced  in  magnitude  (Thornton  2000).  Peters  et  al. 
(2010)  recommend  exercising  caution  when  adjusting  material  properties 
to  expedite  run  times,  especially  when  examining  dynamic  problems 
where  the  fundamental  frequency  of  the  mass  can  be  an  issue. 

At=2]jf  (13> 


where: 

At  =  the  critical  time-step 
m  =  the  particle  mass 
k  =  the  particle  stiffness 

3.3  Horner,  Peters,  and  Carrillo  model 

3.3.1  Model  development 

The  Horner  et.  al  (2001)  model  simulates  the  3-D  movement  of  spherical 
particles  as  they  interact  with  one  another  and  rigid  solid  objects  defined 
by  finite  elements.  The  development  of  the  model  was  driven  by  studies  for 
modeling  large-scale  deformations  in  systems  of  particles.  The  original 
code  was  modified  with  an  emphasis  on  micro-scale  interactions  that 
could  be  run  on  a  standard  PC  in  serial  based  code.  The  current  code  is 
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written  in  FORTRAN  95  and  is  intended  for  the  micro-mechanical 
research  of  engineering  materials. 

The  objective  of  Horner  et.  al  (2001)  was  the  development  of  a  3-D  large- 
scale  soil  modeling  system.  The  implementation  of  a  discrete  element 
method  (DEM)  model  for  the  simulation  of  vehicle  soil  interaction  was  the 
focus  of  the  effort.  The  developed  model  consisted  of  three  sections  to 
model  the  physics  of  soils  to  include  calculation  of  boundary  and  global 
body  forces,  calculate  forces  between  particles,  and  integrate  equations  of 
motion  in  accordance  with  Newton’s  second  law.  The  original  algorithm 
loops  through  the  previously  described  actions  on  a  particle  by  particle 
basis  for  the  duration  of  the  simulation.  Horner’s  model  has  been 
optimized  through  code  efficiencies  and  a  new  contact  detection 
algorithm.  A  serial  version  of  this  optimized  code  is  employed  in  this 
study.  As  of  this  study,  due  to  the  computational  time  requirements,  only  a 
few  tens  of  thousands  of  particles  can  be  simulated  on  a  personal 
computing  platform.  One-to-one  simulation  of  a  triaxial  compression 
specimen  composed  of  sand  would  require  approximately  200,000 
particles.  Simulation  of  field  applications  using  the  same  sand  would 
require  well  over  1  million  particles. 

The  authors  (Horner  et.  al  2001)  investigated  the  use  of  coarsened 
particulate  models  where  each  particle  in  the  simulation  represents  a  large 
volume  of  particulate  media.  DEM  was  used  to  simulate  the  movement  of 
a  fined  plow  through  a  mass  of  sand.  The  authors  determined  that  the 
location  and  width  of  the  shear  band  as  well  as  the  motion  of  particles 
within  the  band  could  be  accurately  predicted  with  the  DEM  model. 

3.3.2  Elements  of  the  DEM  model 

The  DEM  model  used  in  this  study  is  composed  of  spherical  particles  and 
rigid  objects.  Other  elements  that  are  available  in  the  model  include  non- 
spherical  particles  and  membranes.  The  DEM  model  was  developed  to 
provide  data  on  inter-  and  intra-  particle  forces  and  displacements  of  the 
particles.  Therefore,  the  other  objects  are  limited  in  the  manner  in  which 
they  interact.  The  particles  can  interact  with  both  particles  and  other 
objects  in  the  model;  however,  rigid  objects  and  membrane  can  only 
interact  with  particles  and  themselves. 
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3.3.3  Spherical  particles 

Spherical  particles  are  the  primary  elements  of  interest  in  the  current 
study.  Although  non-spherical  particles  more  accurately  represent  the 
shape  of  actual  soils,  spherical  particles  were  selected  to  reduce  the 
computation  time  for  the  model.  The  particles  are  defined  by  the  location 
of  the  centroid,  radius,  and  density  of  the  particle.  The  state  of  the  particle 
at  any  time  increment  is  defined  by  the  position,  linear  and  rotational 
velocity,  and  cumulative  force  and  moment.  The  cumulative  forces  and 
moments  arise  as  each  particle  interacts  with  other  particles  and  rigid 
bodies  in  the  model.  Additionally,  the  mass  of  the  particle  is  assumed  to  be 
uniformly  distributed  so  that  the  moment  of  inertia  tensor  is  isotropic 
(Peters  et  al.  2010). 

3.3.4  Facets 

The  interaction  of  particles  with  rigid  objects  in  the  model  is  defined 
through  contact  with  facets.  Facets  are  3-D,  oriented  triangular  surfaces 
used  to  define  rigid  objects  (Figure  4).  The  facets  are  defined  similarly  to 
finite  elements  using  connectivity  tables  to  monitor  the  coordinates  of 
each  facet’s  vertices.  A  particle  is  considered  to  be  in  contact  with  a  facet  if 
the  distance  between  the  particle  centroid  and  the  facet  face  is  less  than 
the  radius  of  the  particle,  and  the  potential  contact  point  between  the 
particle  and  facet  lies  within  the  boundaries  of  the  facet.  These  contact 
criteria  are  insufficient  at  the  edge  or  corner  of  a  rigid  object  where  the 
potential  contact  point  can  lie  outside  the  boundary  of  the  facet  yet  the 
particle  and  facet  are  in  contact.  Line  elements  are  used  in  the  model  to 
account  for  occurrences  where  the  particle  is  in  contact  with  the  edge  or 
corner  of  an  object. 
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Figure  4.  Illustration  of  a  discrete  element  facet. 


3.3.5  Lines 

The  line  element  is  used  to  connect  two  facets  at  the  edge  or  corner  of  an 
object  (Figure  5).  Line  elements  are  defined  similarly  to  facets  except  the 
connectivity  only  requires  the  two  nodes  defining  the  endpoints.  The 
contact  criteria  for  facet  elements  is  insufficient  at  the  edge  or  corner  of  an 
object  where  a  particle  can  be  in  contact  with  a  facet;  however,  the 
potential  point  of  contact  lies  outside  of  the  facet  boundary.  The  line 
element  is  used  to  resolve  the  inherent  ambiguity. 

Figure  5.  Illustration  of  discrete  element  line  element. 


3.3.6  Rigid  objects 

Rigid  objects  are  the  principal  means  for  imparting  boundary  conditions 
on  a  system  of  particles.  The  rigid  object  is  defined  by  the  collection  of 
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facets  and  lines  that  make  up  the  surface,  the  mass  of  the  object,  and  the 
moment  of  inertia  (Figure  2  and  Figure  6).  Rigid  objects  can  move  through 
three  linear  displacement  and  three  rotational  degrees  of  freedom.  The 
movement  of  the  rigid  object  can  be  controlled  using  either  force  (or 
moment)  or  velocity  (displacement  and  rotation)  variables.  The  movement 
of  the  rigid  objects  is  tracked  by  updating  the  nodes  of  the  composing 
facets  and  lines  using  a  coordinate  transformation  consisting  of  translation 
and  rotation  about  the  center  of  mass  (Peters  et  al.  2010).  The  forces  and 
moments  applied  to  the  facets  are  accumulated  at  the  center  of  mass  for 
the  six  degrees  of  freedom.  Unlike  finite  elements,  the  facet  and  lines  are 
treated  independently  and  therefore  do  not  have  to  be  connected. 

Figure  6.  Illustration  of  a  rigid  object  (plow)  rendered  in  the  discrete  element  model. 


3.4  Parametric  study 

DEM  models  allow  for  simulation  of  particulate-based  materials  in  a  wide 
variety  of  engineering  applications.  The  primary  parameters  that  control 
the  behavior  of  a  system  of  compacted  particles  include:  specific  gravity, 
contact  stiffness,  surface  friction,  and  rolling  resistance.  While  specific 
gravity,  contact  stiffness,  and  surface  friction  are  properties  of  the  mod¬ 
eled  media,  rolling  resistance  is  a  calibration  factor  employed  to  account 
for  the  use  of  spherical  particles  in  modeling  non-spherical  materials. 
When  simulating  triaxial  compression  testing,  the  axial  force,  confining 
pressure,  and  loading  rate  can  be  specified. 

Prior  to  beginning  the  research  effort  outlined  in  Chapter  1,  a  parametric 
study  was  undertaken  to  gain  familiarity  with  the  influence  of  the  primary 
material  parameters  on  engineering  behavior.  The  interaction  of  particles 
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is  modeled  using  sets  of  springs,  dashpots,  and  a  frictional  slider.  An  effort 
was  made  to  accurately  represent  the  materials  of  interest  as  recom¬ 
mended  by  O’Sullivan  et.  al  (2004).  Therefore,  material  properties 
representative  of  a  USCS  classified  (SP)  quartz  sand  were  used  during  the 
parametric  study  to  include  a  specific  gravity  of  2.65,  normal  and  shear 
spring  stiffness  ( Kn  and  Xs )  of  5.0  x  106  psi,  friction  coefficient  (u)  of  0.50, 
and  rolling  resistance  ( R )  of  50.0  to  account  for  the  use  of  spherical 
particles.  Additionally,  the  triaxial  testing  was  conducted  under  strain 
control  with  an  initial  deviator  stress  of  50  psi  and  a  loading  rate  of 
0.05  in. /sec.  The  influence  of  particle  stiffness,  surface  friction,  and  rolling 
resistance  were  investigated  by  varying  each  property  within  a  finite  range 
while  all  other  parameters  were  held  constant. 

3.4.1  Effect  of  particle  stiffness 

The  effect  of  particle  stiffness  on  the  engineering  behavior  of  compacted 
cubical  specimens  subjected  to  triaxial  compression  was  investigated  using 
a  serial  version  of  the  Horner,  Peters,  and  Carrillo  DEM  model.  The  speci¬ 
mens  consisted  of  approximately  7,000  uniform  particles  with  a  diameter 
of  0.80  in.  and  the  properties  of  SP  sand  as  previously  described.  The 
specimen  generation  process  is  discussed  in  detail  in  Chapter  5. 

Particles  interact  in  the  DEM  model  in  accordance  with  linear-elastic  con¬ 
tact  laws  based  upon  Hooke’s  law  rather  than  a  Hertzian  contact  law 
typical  in  most  DEM  soil  mechanics  models.  The  particle  stiffness  is 
determined  by  the  magnitude  of  the  spring  constant  assigned  in  the 
normal  (Xn)  and  shear  (Xs)  directions.  For  this  research,  the  ratio  of 
normal  to  shear  stiffness  was  maintained  at  1:1  and  illustrates  the  stress- 
strain  relationship  on  the  particle  matrix  for  this  stiffness  condition  for  K 
values  ranging  from  500,000  to  50,000,000  psi  (shown  in  Figure  7). 
Cundall  and  Strack  (1979)  observed  the  ratio  of  normal  and  shear  stiffness 
to  have  a  significant  influence  on  the  behavior  of  a  particulate  system.  The 
influence  however,  was  reduced  at  low  friction  angles  due  to  the  early 
onset  of  slip.  Belheine  et  al.  (2009)  also  observed  that  elastic  response  is 
controlled  by  the  ratio  of  shear  to  normal  contact  stiffness.  The  shear  force 
is  controlled  exclusively  by  the  normal  force  once  slip  has  occurred  in 
accordance  with  the  Mohr-Coulomb  criteria  (Cundall  and  Strack  1979). 
Due  to  the  use  of  spherical  particles,  it  was  hypothesized  that  the  onset  of 
slip  would  occur  early  in  testing  and  the  use  of  a  1:1  ratio  was  found  to 
approximate  that  response  effectively  in  the  DEM. 
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The  strength  of  the  particulate  system  modeled  using  DEM  was  greatly 
influenced  by  the  prescribed  contact  stiffness.  Intuitively,  as  the  stiffness  is 
increased  the  peak  stress  also  increases  when  the  strain  is  held  constant. 

In  particular,  an  increase  in  stiffness  results  in  a  significant  increase  in 
shear  stress.  The  mean  stress  measured  during  triaxial  compression  is  also 
increased  but  to  a  lesser  degree.  In  addition,  an  increase  in  stiffness  results 
in  greater  dilative  tendency  as  indicated  by  significant  drops  in  calculated 
pore  water  pressure. 


Figure  7.  Effect  of  particle  stiffness  on  stress-strain  behavior. 


3.4.2  Effect  of  surface  friction 

Similar  to  particle  stiffness,  the  influence  of  surface  friction  on  the  engi¬ 
neering  behavior  of  particulate  assemblies  was  investigated  using  the  DEM 
model.  The  coefficient  of  friction  (u)  describes  the  force  of  friction  between 
two  contacting  particles  in  accordance  with  the  Mohr-Coulomb  criterion 
(u  =  tan  cp).  Particle  stiffness  and  rotational  resistance  were  held  constant 
as  the  particle  surface  friction  coefficient  was  varied  from  o.i  to  5.0  as 
shown  in  Figure  8,  noting  that  during  initial  loading  the  particle 
rearrangement  caused  shear  stress-strain  inconsistencies.  Increasing 
inter-particle  friction  results  in  increased  initial  stiffness  (Cui  et.  al  2007). 
Friction  at  the  contact  controls  the  macro-deformation  and  strength 
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behavior  of  granular  materials  (Oda  and  Iwashita  1999).  Deformation 
depends  strongly  on  local  friction  (Belheine  et  al.  2009).  The  friction 
activated  at  the  contacts  between  structural  units  determines  the  shear 
strength  of  a  particulate  material.  The  mechanical  behavior  is  closely 
related  to  the  translational  and  rotational  friction  (Feda  1982).  However, 
Cavarretta  et.  al  (2010)  found  macro-scale  performance  to  be  more  greatly 
influenced  by  particle  shape  than  surface  roughness. 

Overall  shearing  resistance  is  increased  with  greater  inter-particle  friction. 
The  application  of  pore  pressures  can  be  used  to  overcome  inter-particle 
friction  and  force  a  particulate  material  into  a  fluid  state  (Feda  1982).  At 
high  stress  levels,  the  influence  of  surface  friction  is  reduced  (Cavarretta 
et.  al  2010).  The  magnitude  of  stress  fluctuations  is  influenced  by  system 
size  and  particle  roughness  (Shearer  and  Schaeffer  1997). 


Figure  8.  Effect  of  particle  surface  friction  on  stress-strain  behavior. 


Figures  7  and  8  exhibit  stress-strain  behaviors  of  simulated  specimens 
undergoing  triaxial  compression  testing.  The  peak  shear  stress  of  the 
compacted  system,  and  likewise  the  peak  mean  stress  increased  with 
increased  inter-particle  friction,  u.  However,  the  inter-particle  friction  is 
not  proportional  with  bulk  material  properties.  The  influence  of  inter¬ 
particle  friction  is  greater  at  lower  bulk  values;  at  higher  bulk  values,  an 
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equal  change  in  inter-particle  friction  has  a  reduced  effect.  At  sufficiently 
low  values  of  friction,  the  particle  assembly  was  observed  to  rotate  freely, 
similar  to  a  ball  bearing  and  the  system  reaches  a  steady  state  of 
consolidation. 

3.4.3  Effect  of  rolling  resistance 

The  influence  of  rolling  resistance  was  also  investigated  using  the  DEM 
model.  Rolling  resistance  can  be  applied  when  modeling  irregularly 
shaped  granular  materials  using  spherical  particles.  Resistance  to  rolling  is 
not  a  material  property.  Instead,  it  is  a  calibration  parameter  used  in 
modeling  of  spheres  to  ascertain  the  engineering  response  in  a  pseudo 
non-spherical  manner.  The  rolling  resistance  coefficient  (R)  is  a  scaling 
factor  that  increases  the  rotational  and  torsional  moments  to  account  for 
the  resistance  that  would  be  provided  by  the  geometry  of  irregular  shaped 
particles  according  to  Equation  14. 


m-  =  m-  +(RxAt  x  Aw-  )  (14) 

“Using  rolling  resistance  between  spherical  discrete  elements  is  justified 
since  real  grains  are  generally  not  spherical  and  may  exhibit  a  rough  sur¬ 
face  texture  which  can  be  covered  with  a  thin  film  of  weathered  products.” 
(Belheine  et  al.  2009).  Rolling  resistance  has  great  influence  on  the  con¬ 
tact  behavior  of  the  granular  material.  Rolling  resistance  does  not  affect 
the  translational  motion  but  does  affect  the  angular  motion  (Belheine  et  al. 
2009).  Feda  (1982)  states  that  if  contact  stiffness  is  increased  to  prevent 
sliding  and  rotation,  a  particle  assembly  will  lose  freedom  of  motion  and 
behave  as  a  continuous  solid  which  would  create  an  unrealistic  particulate 
model.  Therefore,  while  particle  stiffness  and  the  surface  friction 
coefficient  were  held  constant,  the  rolling  resistance  factor  was  varied 
from  5.0  to  500.0  as  shown  in  Figure  9. 
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Figure  9.  Effect  of  particle  rolling  resistance  on  stress-strain  behavior. 
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The  influence  of  rolling  resistance  was  examined  using  a  simulated  sand 
material  undergoing  triaxial  compression  testing.  Increasing  rolling  resis¬ 
tance  results  in  a  minor  increase  in  peak  stress.  In  particular,  the  shear 
stress  in  the  system  is  slightly  increased  with  essentially  no  increase  in 
mean  stress.  However,  the  greatest  influence  of  increasing  rolling  resis¬ 
tance  can  be  seen  in  the  volumetric  changes  induced  in  the  material 
undergoing  shearing.  Increasing  rolling  resistance  increases  contractency 
and  decreases  dilatancy  potential.  In  Figure  9,  at  a  rolling  resistance  of  5, 
the  sand  material  exhibits  a  lower  overall  shear  strength  than  at  a  value  of 
500.  This  increase  in  strength  at  a  high  rolling  resistance  supports  the 
occurrence  of  increased  dilatancy  between  particle  contacts. 
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4  Experimental  Validations 

4.1  Testing  program 

The  objectives  of  the  physical  testing  regimen  were  to  develop  calibration 
curves  for  the  DEM  and  to  validate  the  mechanical  performance  of 
bimodal  particle  distributions  subjected  to  simulated  triaxial  compression 
testing.  Laboratory  testing  consisted  of  triaxial  compression  of 
consolidated  particulate  specimens  in  accordance  with  ASTM  D  4767, 
Standard  test  method  for  consolidated  undrained  triaxial  compression 
test  for  cohesive  soils  (ASTM  2004).  The  standard  testing  method  is  often 
used  to  determine  the  shear  strength  of  soils  in  a  variety  of  engineering 
applications  to  include  embankment  stability  analysis,  earth  pressure 
calculations,  and  foundation  design  (ASTM  2004).  Undrained  testing  was 
selected  for  the  granular  medium  as  pore  pressure  response  allows 
measurement  of  volume  response  in  a  more  comprehensive  and  finer 
resolution  than  vertical  and  horizontal  strain  measurements  of  the 
cylindrical  sample.  An  effective  stress  path  can  be  calculated  from  the  CD 
DEM  model  to  allow  comparisons  to  the  CU  laboratory  results.  Initially, 
assemblies  of  uniform  particle  sizes  were  constructed  and  tested  to 
produce  data  for  calibration  of  the  DEM  model.  Bimodal  blends  were  then 
used  to  validate  the  behavior  observed  in  the  numerical  simulations.  Both 
stainless  steel  and  polypropylene  materials  were  used  with  sphere 
diameters  of  0.75  in.  and  0.1875  in.  Peters  and  Berney  (2010)  used  an 
identical  testing  approach  to  examine  the  behavior  of  bimodal  blends  of 
sand  and  silty-clay.  The  central  theme  of  that  work  was  the  identification 
of  a  critical  mixture  ratio  identified  using  percolation  theory. 

4.2  Idealized  particles 

The  particles  used  in  the  laboratory  investigation  consisted  of  0.75-in.-  and 
o.i875-in.-diam  spheres  composed  of  stainless  steel  or  polypropylene  as 
shown  in  Figure  10.  All  the  spheres  used  in  this  study  were  obtained  from 
Bal-tec.  Bal-tec  is  a  division  of  Micro  Surface  Engineering  Incorporated, 
located  in  Los  Angeles,  California,  specializing  in  the  manufacture  of 
precision  balls  and  bearings.  The  stainless  steel  balls  are  type  440C  stainless 
steel  that  consists  primarily  of  chromium  and  ferrite.  Type  440C  stainless 
steel  is  used  widely  in  anti-friction  applications  where  high  hardness  and 
corrosion  resistance  are  required.  The  material  has  Rockwell  Hardness  of 
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58  and  a  density  of  0.277  lb/in.3.  Polypropylene  is  a  thermoplastic  polymer 
used  in  a  number  of  engineering  and  industrial  applications.  The  material  is 
extremely  corrosion  resistant  but  is  flexible  and  resists  fatigue.  The 
polypropylene  material  used  in  the  manufacture  of  the  spheres  has  an 
average  Rockwell  Hardness  of  90  and  a  density  of  0.033  lb/in.3.  The 
manufacturing  quality  of  the  spheres  is  indicated  by  the  respective  grades. 
The  stainless  steel  spheres  were  a  grade  25,  indicating  maximum  deviation 
from  target  diameter  and  sphericity  of  2.5  x  io_5  in.  The  polypropylene 
spheres  were  a  grade  3  with  maximum  deviation  of  3  x  icr6  in. 


Figure  10.  Particles  examined  in  study:  stainless  steel  (left) 
and  polypropylene  (right). 


4.2.1  Consolidated-undrained  triaxial  testing 

Consolidated-undrained  triaxial  compression  testing  conducted  in  accord¬ 
ance  with  ASTM  D  4767  (ASTM  2004)  provides  data  on  the  strength  and 
stress-strain  relationship  of  a  cylindrical  particulate  specimen.  Specimens 
were  saturated,  consolidated  isotropically,  and  sheared  using  a  constant 
rate  of  compression  without  drainage.  Determination  of  total  and  effective 
stresses  is  possible  through  measurement  of  axial  load,  axial  deformation, 
and  pore-water  pressure. 

4.2. 1.1  Testing  apparatus 

The  testing  apparatus  used  in  this  study  was  custom  manufactured  by 
MTS  Systems  Corporation  for  the  U.S.  Army  Engineer  Research  and 
Development  Center  (ERDC).  The  apparatus  consisted  of  an  axial  loading 
frame,  hydraulic  supply,  electronic  controller,  load  cell,  linear  variable 
differential  transformer  (LVDT)  and  pore-water  pressure  transducer,  all 
conforming  to  the  standards  provided  in  ASTM  D  4767  (ASTM  2004)  and 
presented  in  Figure  11. 
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Figure  11.  MTS  triaxial  compression  testing  apparatus. 
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The  triaxial  compression  chamber  was  modified  from  its  initial 
manufacture  to  accommodate  the  instrumentation  required  for  the 
experiments.  The  hydraulically  driven  loading  frame  consists  of  a  MTS 
Model  312  load  frame  with  55-kip  maximum  capacity,  Model  510 
hydraulic  supply  and  Model  293  hydraulic  manifold  capable  of  supplying 
3,000  psi  of  pressure,  and  Flex  Test  GT  digital  controller.  The  frame  can 
be  either  load  or  displacement  controlled.  Displacement  control  was  used 
with  a  loading  rate  of  0.05  in. /min  for  this  investigation.  The  BLH 
Electronics  Model  U3G1  electronic  load  cell  has  a  maximum  capacity  of 
5,000  lbf  and  provides  1  x  io~3  lbf  precision.  The  RDP  Group  Model 
LDC1000A  LVDT  measures  the  axial  deformation  of  the  specimen  with  a 
maximum  vertical  displacement  of  approximately  3.0  in.  and  provides  1  x 
io-7-in.  precision.  The  pore-water  pressure  is  measured  using  a  Sensotec 
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K2  electronic  pressure  transducer  with  an  optimal  operating  range  of  5  to 
350  psi,  providing  1  x  io_6-psi  precision.  The  triaxial  compression  chamber 
is  a  sealed  reinforced  chamber  capable  of  withstanding  pressures  up  to 
250  psi.  The  cap  and  base  are  made  of  aluminum,  the  mounting  pedestal  is 
made  of  brass,  and  the  cylindrical  walls  are  made  of  reinforced  Lucite.  The 
base  and  pedestal  provide  ports  for  the  introduction  and  removal  of  gases 
and  liquids  to  the  chamber,  specimen,  or  both. 

4.2. 1.2  Specimen  preparation 

The  specimens  constructed  for  this  investigation  were  cylindrical  with  a 
4-in.-diam  and  a  target  8-in.  height.  Due  to  the  irregular  packing  of  the 
spheres,  some  deviation  in  the  specimen  height  (±0.25  in.)  was  observed. 
Specimen  preparation  can  have  a  dramatic  effect  on  laboratory  perfor¬ 
mance  (Kuo  and  Frost  1997).  To  reduce  potential  for  introduction  of  error 
into  the  testing  program,  the  procedure  provided  in  ASTM  D  4767  (ASTM 
2004)  was  adhered  to  for  the  construction  of  each  tested  specimen.  The 
average  volumetric  properties  for  the  stainless  steel  and  polypropylene/ 
stainless  steel  blends  prepared  for  this  investigation  are  presented  in 
Tables  2  and  Table  3,  respectively. 

The  number  of  particles  required  to  fill  the  4-in.  by  8-in.  target  specimen 
where  calculated  based  upon  volume,  blended  together  (only  for  bimodal 
blends),  and  compacted  in  a  split  cylindrical  mold  in  six  layers  using  a 
tamping  rod.  The  split  cylinder  was  lined  with  a  0.025-in. -thick  neoprene 
membrane  prior  to  addition  of  the  spheres.  Once  the  target  height  was 
achieved,  a  vacuum  was  applied  to  the  specimen,  now  encased  in  the  neo¬ 
prene  membrane,  and  the  split  cylinder  mold  removed.  The  specimen 
preparation  procedure  is  presented  pictorially  in  Figure  12. 


Table  2.  Average  volumetric  properties  of  stainless  steel  blends. 


Particle 

Fraction  % 
(coarse  to  fine) 

15  psi  eff.  confining  stress 

50  psi  eff.  confining  stress 

Mass 

(gm) 

Height 

(in.) 

Void  Ratio 

Mass 

(gm) 

Height 

(in.) 

Void  Ratio 

0-100 

8,300 

7.93 

0.60 

8,300 

7.95 

0.60 

40-60 

9,060 

7.82 

0.43 

9,060 

7.90 

0.43 

45-55 

9,230 

7.98 

0.43 

9,231 

8.04 

0.43 

48-52 

9,174 

7.90 

0.43 

9,174 

8.00 

0.43 

50-50 

9,228 

7.88 

0.42 

9,228 

7.99 

0.42 

52-48 

9,281 

7.87 

0.41 

9,281 

7.99 

0.41 

60-40 

9,244 

7.78 

0.40 

9,244 

7.88 

0.40 
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Table  3.  Average  volumetric  properties  of  polypropylene/ 
stainless  steel  blends. 


Particle 

Fraction  % 
(coarse  to  fine) 

15  psi  eff.  confining  stress 

50  psi  eff.  confining  stress 

Mass 

(gm) 

Height 

(in.) 

Void  Ratio 

Mass 

(gm) 

Height 

(in.) 

Void  Ratio 

0-100 

8,300 

7.93 

0.60 

8,300 

7.95 

0.60 

40-60 

5,888 

8.16 

0.48 

5,888 

8.14 

0.47 

45-55 

5,505 

8.20 

0.45 

5,505 

8.18 

0.45 

48-52 

5,421 

8.08 

0.45 

5,520 

8.10 

0.45 

50-50 

5,149 

8.16 

0.45 

5,149 

8.13 

0.45 

52-48 

4,977 

8.19 

0.44 

4,977 

8.13 

0.44 

60-40 

4,260 

7.98 

0.42 

4,260 

7.97 

0.42 

Figure  12.  Specimen  preparation  procedure:  (a)  lining  mold  with 
membrane,  (b)  prepared  mold,  (c)  removing  split  mold, 

(d)  specimen  ready  for  testing. 


4.2. 1.3  Testing  procedure 

Aluminum  top  and  bottom  caps  were  applied  to  the  specimen  and  the 
assembly  mounted  on  the  pedestal  within  the  triaxial  compression  cham¬ 
ber.  Water  and  vacuum  lines  were  connected  to  allow  water  to  enter  and 
exit  the  specimen  and  a  vacuum  to  be  applied  to  both  the  specimen  and 


ERDC/GSL  TR-16-29 


41 


the  chamber.  The  chamber  was  then  sealed  and  mounted  within  the  load¬ 
ing  apparatus.  Air  was  used  to  apply  the  confining  pressure  in  this  investi¬ 
gation.  The  specimen  was  then  saturated  by  allowing  water  to  flow  from 
the  bottom  cap,  through  the  specimen,  to  the  top  cap.  Once  the  specimen 
was  nearly  saturated,  70  psi  of  backpressure  was  applied  to  force  any 
remaining  air  into  solution  and  fill  all  voids  with  de-aired  water. 

A  check  of  the  B-value,  the  ratio  of  the  change  in  pore-water  pressure 
resulting  from  changes  in  confinement  pressure,  was  then  conducted  to 
determine  the  responsiveness  of  the  specimen.  A  B-value  greater  than  or 
equal  to  0.95  indicated  the  ability  to  begin  testing;  a  value  less  than  0.95 
indicated  that  additional  back  pressurization  was  required.  Holding  the 
backpressure  constant  and  increasing  the  chamber  pressure  was  used  to 
achieve  the  target  effective  confining  pressure.  Confining  pressures  of 
15  psi  and  50  psi  were  used  in  this  investigation.  The  effective  confining 
pressure  was  held  constant,  the  drainage  valves  closed,  and  an  axial  load 
was  applied  to  the  specimen  to  induce  shear.  The  axial  loading  was 
displacement  controlled  at  a  rate  of  0.05  in./min.  Testing  was 
discontinued  when  10  percent  axial  strain  or  shear  failure  of  the  specimen 
was  observed. 

4.2.2  Testing  results 

Data  were  collected  during  triaxial  compression  testing  at  a  sampling  rate 
of  2.0  Hz  and  included  time  elapsed,  axial  load,  axial  displacement,  cham¬ 
ber  pressure,  and  backpressure.  The  data  were  used  to  calculate  shear 
stress,  effective  mean  stress,  shear  strain,  and  excess  pore-water  pressure 
using  Equations  15  through  19,  respectively. 

0  shear  1-°3  (15) 


where: 

0  shear  =  shear  stress 

=  vertical  principal  stress  (axial  load/specimen  area) 

03  =  chamber  pressure 

The  change  in  cross-sectional  area  during  shear  was  approximated  using 
Equation  16. 
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A 


adjusted 


n  d2 


Ml 


(16) 


where: 


Aadjusted 


d 


Saxial 


cross-sectional  area  of  the  specimen  resulting  from 
increasing  axial  strain 
specimen  diameter 

axial  strain  (change  in  height/initial  height) 


8 


shear 


1.5*8 


axial 


(17) 


where  Sshear  is  the  shear  strain. 


A|i  =  n-n,. 


(18) 


where: 


a.|u  =  excess  pore-water  pressure 
p  =  backpressure  at  any  time 

Hi  =  initial  backpressure 


o„,„„„  = 


o1  +2o^ 


h 


where  0mean  is  the  effective  mean  stress. 


(19) 


The  initial  set  of  experiments  was  conducted  on  uniform  assemblies  of 
o.i875-in.-diam  spheres.  The  testing  was  conducted  on  both  the  stainless 
steel  and  polypropylene  materials  at  15-psi  and  50-psi  confining  pressure. 
The  objective  of  these  experiments  was  to  develop  characteristic  curves  for 
the  materials  used  in  the  study  and  apply  these  curves  for  the  calibration 
of  the  DEM  material  model.  Two  tests  were  run  for  each  material  type  and 
confining  stress  level  combination  and  the  resulting  shear  stress,  excess 
pore  water  pressure  and  effective  mean  stress  for  each  combination  were 
averaged  together  for  each  level  of  shear  strain.  The  average  stress-strain 
characteristic  curves  developed  using  the  procedure  previously  described 
are  presented  in  Figure  13. 
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Figure  13.  Stress-strain  behavior  of  uniform  0.1875-in.  spheres. 
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Once  the  characteristic  curves  were  developed,  blends  of  the  stainless  steel 
spheres  were  created  and  tested  to  examine  the  influence  of  particle-size 
distribution  on  engineering  performance.  The  specimens  were  generated 
and  tested  as  described  in  the  “Consolidated-undrained  triaxial  testing” 
section.  Initially,  blends  of  40  percent  (by  volume)  o.75-in.-diam  spheres 
and  60  percent  o.i875-in.-diam  spheres  were  examined  at  both  15-psi  and 
50-psi  confining  pressure.  These  blends  were  designated  40  to  60  percent 
(0.75  to  0.1875  in.).  Two  tests  were  run  at  each  stress  level  and  confining 
stress  level  combination  and  the  resulting  shear  stress,  excess  pore  water 
pressure  and  effective  mean  stress  for  each  combination  were  averaged 
together  for  each  level  of  shear  strain. 

This  process  was  repeated  for  increasing  proportions  of  0.75-in.  (coarse) 
particles.  The  designations  for  the  examined  particle-size  distributions 
include  45-55, 48-52,  50-50,  52-48,  and  60-40  (coarse  percentage-fine 
percentage)  as  presented  in  Table  4.  The  particle-size  distributions  used  in 
this  study  were  selected  to  identify  the  proportions  at  which  the  specimen 
behavior  is  dominated  by  one  proportion  or  the  other  according  to 
percolation  theory  (Peters  and  Berney  2010).  The  stress-strain  plots  for 
the  triaxial  tests  run  in  this  investigation  at  15  psi  and  50  psi  are  presented 
in  Figures  14  and  15,  respectively. 
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O’Sullivan  et.  al  (2006)  observed  a  stiff  initial  response  until  the  peak 
strength  was  achieved,  followed  by  localization  formation  at  the 
origination  of  shear  planes  and  dilation  in  triaxial  testing  of  ideal  spheres. 
Similar  results  were  obtained  in  this  study  as  shown  in  Figure  14  through 
Figure  17.  The  stainless  steel  assemblies  exhibited  stiff  initial  response 
followed  by  localization  and  dilation.  In  this  study,  at  a  confining  stress 
level  of  15  psi  the  stiffest  blends  were  the  40-60,  the  weakest  were  the 
60-40,  and  the  blends  approaching  the  percolation  threshold  purposed  by 
Peters  and  Berney  (2010)  fell  in  the  middle.  At  50-psi  confining  pressure 
the  blends  40-60  and  60-40  exhibited  the  least  shear  strength  and  the 
blends  approaching  the  percolation  threshold  displayed  greatest  stiffness. 
Slight  differences  in  behavior  with  increasing  proportions  of  large  particles 
were  observed. 


Table  4.  Target  number  of  spheres  for  consolidated  laboratory  specimen. 


Particle  Fraction  % 
(coarse  to  fine) 

Number  of  0.75-in.-diam 
Particles 

Number  of  0.1875-in.-diam 
Particles 

0-100 

0 

19,091 

40-60 

126 

12,500 

45-55 

148 

11,475 

48-52 

153 

11,023 

50-50 

162 

10,568 

52-48 

171 

10,114 

60-40 

198 

8,295 

At  the  conclusion  of  testing  on  bimodal  blends  of  stainless  steel  spheres, 
the  testing  protocol  was  extended  to  bimodal  blends  of  polypropylene  and 
stainless  steel  sphere  mixtures.  The  large  proportion  (0.75-in.  diam)  of 
these  blends  was  composed  of  polypropylene  spheres  and  the  small  pro¬ 
portion  (0.1875-in.  diam)  was  stainless  steel  spheres.  The  previous  series 
of  tests  focused  on  the  influence  of  particle-size  distribution,  while  the 
polypropylene/stainless  steel  series  attempted  to  examine  the  influence  of 
particle  stiffness.  Many  common  engineering  materials  including  soils  and 
asphalt  concrete  are  composed  of  materials  with  considerably  different 
stiffness  properties.  An  identical  testing  procedure  was  conducted  for  the 
polypropylene/stainless  steel  blends  as  described  for  the  stainless  steel 
bimodal  blends.  The  stress-strain  plots  for  the  triaxial  tests  run  in  this 
investigation  at  15  psi  and  50  psi  are  presented  in  Figure  16  and  Figure  17, 
respectively  as  averages  of  two  replicates  conducted  for  each  mixture. 
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The  bimodal  blends  of  polypropylene  and  stainless  steel  spheres  exhibited 
similar  performance  trends  as  compared  to  the  assemblies  of  only 
stainless  steel  spheres  with  stiff  initial  response,  followed  by  localization 
formation,  and  dilation.  However,  due  to  the  reduced  stiffness  of  the  large 
proportion  composed  of  polypropylene  spheres,  increasing  the  large  pro¬ 
portion  resulted  in  a  decrease  of  both  strength  and  stiffness.  Most  of  the 
specimens  reached  10  percent  axial  strain  without  failure.  The  behavior 
exhibited  by  the  polypropylene/stainless  steel  blends  followed  expectation. 
The  blends  with  a  greater  proportion  of  the  large  fraction  exhibited 
reduced  stiffness  and  a  reduced  propensity  to  dilate  with  the  60-40  blend 
exhibiting  continued  consolidation  through  10  percent  shear  strain,  then 
switching  to  a  dilative  state.  Also  dissimilar  from  the  stainless  steel  only 
tests,  a  distinct  difference  in  behavior  can  be  observed  with  increasing 
large  (soft)  proportion. 


Figure  14.  Stress-strain  behavior  of  stainless  steel  spheres  at  15  psi. 


- 0-100 

- 40-60 

- 45-55 

- 48-52 

- 50-50 

- 52-48 

- 60-40 


ERDC/GSL  TR-16-29 


46 


Figure  15.  Stress-strain  behavior  of  stainless  steel  spheres  at  50  psi. 
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Figure  16.  Stress-strain  behavior  of  polypropylene/ 
stainless  steel  spheres  at  15  psi. 
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Figure  17.  Stress-strain  behavior  of  polypropylene/ 
stainless  steel  spheres  at  50  psi. 
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4.3  Naturally  occurring  particles 

The  use  of  uniform  spheres  allows  for  the  generation  of  numerical  models 
that  accurately  represent  the  simulated  particle  geometry.  Numerical  sim¬ 
ulations  provide  valuable  data  on  the  micromechanics  of  particulate  sys¬ 
tems  that  is  unattainable  using  standard  laboratory  measures.  However, 
naturally  occurring  particles  that  are  commonly  used  in  engineering  appli¬ 
cations  deviate  greatly  from  easily  represented  ideal  shapes  (spheres, 
cubes  etc.).  The  use  of  ideally  shaped  particles  (spheres)  to  investigate  the 
micromechanical  response  of  naturally  occurring  particles  was  investi¬ 
gated  in  this  study.  Although  the  geometry  of  the  particles  is  not  accurately 
represented,  the  application  of  the  DEM  model  can  still  provide  valuable 
information  on  the  behavior  of  particulate  assemblies  of  interest.  Labora¬ 
tory  testing  of  naturally  occurring  materials  was  not  pursued  as  part  of  this 
investigation.  However,  laboratory  data  published  by  Peters  and  Berney 
(2010)  on  bimodal  particle  distributions  of  sand  and  silty-clay  were 
examined  and  provided  a  platform  to  calibrate  the  DEM  model  and 
validate  the  laboratory  observed  engineering  behavior.  This  process 
provided  the  necessary  guidance  to  calibrate  the  polypropylene  and  steel 
materials  used  in  this  study. 
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Peters  and  Berney  (2010)  investigated  the  constitutive  response  of  gap- 
graded  binary  mixtures  at  critical  mixture  ratios.  Important  engineering 
materials  including  flexible  pavements,  low  permeability  retaining  struc¬ 
tures,  and  filters  are  often  approximated  as  gap-graded  particle-size  dis¬ 
tributions.  The  study  focused  on  a  laboratory  study  of  gap-graded 
materials  analyzed  using  modern  micro-mechanical  research  techniques. 
The  central  theme  of  the  work  was  the  identification  of  a  critical  mixture 
ratio  identified  using  percolation  theory. 

4.3.1  Summary  of  testing  presented  by  Peters  and  Berney 

The  mechanical  behavior  of  binary  mixtures  composed  of  sand  and  silty- 
clay  as  the  relative  fraction  of  each  material  was  examined  by  Peters  and 
Berney  (2010).  Triaxial  compression  testing  was  conducted  on  1.4-in.- 
diam  by  3.o-in.-tall  cylindrical  specimens  of  USCS  classified  SP  sand  and 
CL  silty-clay.  The  specimens  were  compacted  using  a  pneumatic  kneading 
compactor  with  tip  forces  ranging  from  5  to  60  psi.  The  specimens  were 
tested  in  a  triaxial  apparatus  in  accordance  with  ASTM  standard  test 
method  D  4767  (ASTM  2004).  The  researchers  found  evidence  of  the 
percolation  threshold  at  sand  fractions  between  45  and  48  percent. 
Mixtures  having  reached  the  threshold  displayed  erratic,  even  unstable 
behavior  with  respect  to  expected  mechanical  response.  The  traditional 
relationship  between  mean  stress  and  void  ratio  based  on  the  critical  state 
theory  was  also  observed  to  break  down  as  the  coarse  material  nears  its 
threshold  fraction  (Peters  and  Berney  2010).  Additionally,  the  sand 
fraction  controlled  the  achievable  compaction  in  the  silty-clay  fraction. 
The  authors  concluded  that  the  observed  behavior  was  due  to  the 
development  and  evolution  of  force  chains  and  not  a  reduction  in  void 
ratio  caused  by  increased  confining  pressure. 

4.3.2  Testing  results 

Results  of  laboratory  testing  as  presented  by  Peters  and  Berney  (2010) 
were  used  to  calibrate  and  validate  DEM  simulations  of  assemblies  of  par¬ 
ticles  representing  a  bimodal  sand/silty-clay  blend.  The  stress-strain 
response  of  binary  blends  of  sand  and  silty-clay  at  15  psi  and  50  psi  was 
reproduced  using  data  from  Peters  and  Berney  (2010)  and  is  presented  in 
Figures  18  and  19.  Characteristic  curves  developed  using  uniform  assem¬ 
blies  of  sand  and  silty-clay  were  used  to  calibrate  the  DEM  model  prior  to 
the  simulation  of  sand/silty-clay  bimodal  blends  and  are  presented  in 
Figure  18  through  Figure  20. 
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Figure  18.  Stress-strain  behavior  of  sand  and  silty-clay  blends  at 
15  psi  used  in  Peters  and  Berney  (2010)  study. 


Figure  19.  Stress-strain  behavior  of  sand  and  silty-clay  blends  at 
50  psi  used  in  Peters  and  Berney  (2010)  study. 
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Figure  20.  Stress-strain  behavior  of  sand  and 
silty-clay  calibration  curves. 
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5  Numerical  Simulation  of  Bimodal  Systems 

5.1  Testing  program 

The  objective  of  the  numerical  simulations  was  to  develop  assemblies  of 
consolidated  spheres  that  exhibit  similar  mechanical  behavior  as  the  spe¬ 
cimens  examined  in  laboratory  testing  and  to  use  the  specimens  to  inves¬ 
tigate  the  stress  transfer  mechanisms  occurring  in  bimodal  particulate 
systems  undergoing  shear.  The  intention  was  not  to  develop  the  most 
physically  realistic  numerical  model  but  to  build  one  from  which  meaning¬ 
ful  data  on  the  particulate  mass  could  be  obtained.  Simulations  were  made 
of  varying  particle-size  distributions  composed  of  stainless  steel  and  poly¬ 
propylene/stainless  steel  spheres  as  outlined  in  Chapter  4.  Simulation  of 
natural  soil  materials  composed  of  bimodal  blends  of  sand  and  silty-clay 
was  also  attempted  using  data  obtained  from  Peters  and  Berney  (2010). 

The  DEM  model  developed  by  Horner  et.  al  (2001)  was  used  to  generate 
particulate  specimens  and  simulate  the  consolidated-undrained  triaxial 
compression  testing  of  cubical  assemblies  of  particles.  The  DEM  model 
consisted  of  a  finite  modeling  space,  confining  walls,  and  approximately 
7,000  spherical  particles.  Using  any  more  particles  in  the  DEM  model 
proved  too  time  consuming  and  difficult  for  the  computer  systems  to 
process  and  so  a  compromise  was  made  at  7,000  to  enable  this  initial 
study  to  continue.  The  finite  modeling  space  allowed  for  tracking  of  the 
model  elements  between  computational  time  increments.  The  confining 
walls  were  composed  of  3-D  facets  defined  similarly  to  finite  elements  and 
representing  rigid  objects.  The  spherical  particles  were  elastic, 
indestructible  elements  defined  by  the  location  of  the  centroid,  particle 
radius,  and  specific  gravity.  The  spherical  particles  are  the  primary 
elements  of  interest  that  provide  data  on  inter  and  intra-particle  forces 
occurring  in  the  consolidated  granular  assembly. 

5.1.1  Specimen  generation 

Consolidated  assemblies  of  spheres  were  subjected  to  simulated  triaxial 
compression  testing  to  determine  the  micro-level  performance  of  varying 
size  distributions  of  particulate  materials.  The  specimens  were  generated 
using  a  growth  algorithm  followed  by  isotropic  consolidation.  Once  initi¬ 
ated,  the  DEM  model  generated  an  assemblage  of  uniform  particles  within 


ERDC/GSL  TR-16-29 


52 


the  model  space  with  diameters  equal  to  the  largest  target  particle  size.  A 
growth  algorithm  then  randomly  selects  particles  in  the  assembly,  based 
upon  the  number  of  particles  required  to  reach  the  desired  size 
distribution,  and  decreases  the  size  of  the  selected  particles  with  each 
computational  time-step  until  the  target  size  is  reached.  Cui  and 
O’Sullivan  (2006)  observed  that  randomly  generating  specimens  in  the 
numerical  environment  results  in  slight  perturbations  in  material  fabric. 
This  effect  can  be  exaggerated  by  the  use  of  low  particle  counts.  However, 
the  researchers  found  that  macro  performance  was  not  significantly 
influence  by  the  minor  differences  in  fabric. 

Concurrent  with  the  change  in  particle  size,  the  six  confining  walls  begin 
moving  towards  the  center  of  the  model  space  under  fixed  velocity  control, 
As  the  rigid  walls  impact  the  particles  contained  within,  the  particles  are 
forced  together  forming  the  consolidated  specimen.  Maximum  density  is 
reached  when  no  changes  in  specimen  dimensions  occur  with  additional 
computational  time-steps. 

The  particle  sizes  used  in  the  DEM  model  were  0.80  and  0.20  in.  diam 
with  size  distributions  similar  to  those  described  in  Chapter  4  and  shown 
in  Table  5.  The  number  and  size  of  particles  used  in  the  numerical  model 
compare  favorably  to  the  ideal  spheres  used  in  physical  testing  (0.75  and 
0.1875  in.).  However,  the  sand  and  silty-clay  grains  used  by  Peters  and 
Berney  (2010)  had  average  diameters  of  0.086  and  0.0001  in.,  respec¬ 
tively.  Their  3.0-  by  1.4-in.  cylindrical  triaxial  specimens  contained  an 
average  4,500  sand  or  7.3  xio10  silty-clay  particles.  The  simulation  of  a 
representative  specimen  of  the  silty-clay  particles  was  not  possible  due  to 
the  prohibitive  computational  expense  (Bathurst  and  Rothenburg  1988). 
Additionally,  one  to  one  correspondence  between  computational  and 
physical  particle  size  is  not  necessary,  because  the  employed  test  bed 
represents  a  representative  elemental  volume  of  particles  that 
approximates  the  mechanical  response  of  the  prototype  medium  rather 
than  being  a  faithful  rendition  of  micromechanical  behavior  (Horner  et.  al 
2001). 


Table  5.  Sphere  distribution  for  numerical  model  consolidated  specimen. 


Particle  Fraction  % 
(coarse  to  fine) 

Number  of  0.80-in.-diam 

Particles 

Number  of  0.20-in.-diam 

Particles 

0-100 

0 

6,859 

40-60 

37 

6,060 

48-52 

72 

6,621 
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Particle  Fraction  % 
(coarse  to  fine) 

Number  of  0.80-in.-diam 

Particles 

Number  of  0.20-in.-diam 

Particles 

50-50 

93 

6,673 

52-48 

106 

6,675 

60-40 

165 

6,672 

5.1.2  Consolidation  of  particles 


Once  a  specimen  of  consolidated  spheres  at  the  maximum  possible  density 
was  generated,  the  next  step  was  to  achieve  the  desired  volumetric  and 
confining  stress  states  prior  to  simulating  triaxial  compression  testing. 
During  the  initial  growth  and  consolidation,  the  spheres  within  the  speci¬ 
men  were  assigned  the  material  properties  of  stainless  steel.  However,  for 
the  polypropylene/stainless  steel  and  sand/silty-clay  blends,  the  material 
properties  of  the  spheres  were  modified  prior  to  transforming  the  speci¬ 
men  further  using  the  values  presented  in  Table  6.  After  the  desired  parti¬ 
cle  properties  were  assigned,  the  confining  walls  were  relaxed  until  the 
desired  confining  stress  was  obtained.  The  walls  were  relaxed  by  switching 
from  velocity  to  force  control  and  specifying  the  force  required  to  reach  a 
confining  stress  of  either  15  or  50  psi.  Evans  et  al.  (2009)  observed  that  it 
is  difficult  to  generate  specimens  with  exact  values  of  both  confining  stress 
and  void  ratio.  The  same  phenomenon  was  experienced  in  this  study.  The 
objective  of  this  investigation  was  to  examine  the  stress  transfer  mecha¬ 
nisms  occurring  in  consolidated  particulate  specimens  that  are  strongly 
influenced  by  the  initial  stress  state.  Therefore,  the  relaxation  process  was 
discontinued  when  the  average  stress  on  the  consolidated  assembly 
reached  the  target  level  ±10  percent. 


Table  6.  Particle  properties  for  materials  simulated  in  DEM  model. 


Stainless  Steel 

Polypropylene 

Sand 

Silty-Clay 

Specific  gravity 

7.80 

0.90 

2.65 

2.70 

Normal  stiffness  (ksi) 

29,000 

290 

6,600 

30 

Shear  stiffness  (ksi) 

29,000 

290 

6,600 

30 

Rolling  resistance 

0.0 

0.0 

7.0 

0.0 

Torsion  resistance 

0.0 

0.0 

7.0 

0.0 

Translational  friction 

0.60 

0.80 

0.70 

1.20 

Rotational  friction 

0.60 

0.80 

0.70 

1.20 

Cui  and  O’Sullivan  (2006)  found  in  constructing  numerical  model  speci¬ 
mens  that  obtaining  an  exact  replica  of  physical  specimens  was  not  possi¬ 
ble.  Typically,  the  void  ratio  of  numerical  specimens  was  slightly  lower 
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than  specimens  constructed  in  the  laboratory.  The  observation  of  Cui  and 
O’Sullivan  (2006)  was  replicated  for  the  coarse  polypropylene/stainless 
steel  and  all  of  the  sand/silty-clay  blends,  but  the  opposite  was  observed 
for  the  stainless  steel  and  fine  polypropylene/stainless  steel  blends  exam¬ 
ined  in  this  study.  Evans  et  al.  (2009)  comment  that  the  relative  density 
across  specimens  should  be  approximately  the  same  even  if  the  void  ratios 
vary  slightly  when  the  specimen  are  generated  and  consolidated  using  an 
identical  method.  An  example  of  a  consolidated  cubical  specimen  is  pre¬ 
sented  in  Figure  21.  The  volumetric  properties  of  the  stainless  steel,  poly¬ 
propylene/stainless  steel,  and  sand/silty-clay  numerical  model  specimens 
are  presented  in  Tables  7  through  Table  9  respectively.  Because  each 
mixture  consisted  of  the  same  number  of  particles,  uniform  mixtures  of  all 
one  size  fraction  (0-100)  exhibited  a  taller  structure,  because  the 
developing  DEM  with  all  one  particle  size  occupies  a  greater  solid  volume 
than  binary  mixtures,  where  many  smaller  particles  were  present  to  reduce 
the  overall  solid  volume. 


Figure  21.  Consolidated  specimen  of  uniform  spheres. 


Table  7.  Volumetric  properties  of  numerical  model  stainless  steel  blends. 


Particle 
Fraction  % 
(coarse  to 
fine) 

15  psi 

50  psi 

Width  (in.) 

Height  (in.) 

Void  Ratio 

Width  (in.) 

Height  (in.) 

Void  Ratio 

0-100 

14.4 

14.5 

0.65 

14.3 

14.4 

0.61 

40-60 

3.6 

3.7 

0.74 

3.3 

3.3 

0.52 

48-52 

3.9 

3.9 

0.68 

3.7 

3.7 

0.54 

50-50 

3.9 

3.9 

0.55 

3.8 

3.8 

0.45 

52-48 

4.2 

4.2 

0.61 

3.9 

3.9 

0.42 

60-40 

4.3 

4.3 

0.44 

4.3 

4.3 

0.43 

ERDC/GSL  TR-16-29 


55 


5.1.3  Simulation  of  triaxial  testing 

The  consolidated-undrained  triaxial  compression  testing  was  simulated 
using  the  DEM  model.  A  zero  volume  change  confining  cell  was  used  to 
approximate  the  effect  of  the  undrained  triaxial  compression  testing.  The 
confining  cell  was  composed  of  six  rigid  bounding  walls  composed  of  facet 
elements.  The  confining  walls  are  frictionless  to  ensure  that  normal 
boundary  induced  stresses  (Belheine  et  al.  2009)  occur  within  the  particu¬ 
late  assembly. 


Table  8.  Volumetric  properties  of  numerical  model  polypropylene/ 
stainless  steel  blends. 


Particle 
Fraction  % 
(coarse  to 
fine) 

15  psi 

50  psi 

Width  (in.) 

Height  (in.) 

Void  Ratio 

Width  (in.) 

Height  (in.) 

Void  Ratio 

0-100 

14.4 

14.5 

0.65 

14.3 

14.4 

0.61 

40-60 

3.6 

3.7 

0.74 

3.3 

3.4 

0.54 

48-52 

3.8 

3.8 

0.60 

3.6 

3.7 

0.50 

50-50 

3.9 

3.9 

0.51 

3.8 

3.8 

0.44 

52-48 

3.9 

3.8 

0.40 

3.8 

3.8 

0.40 

60-40 

4.2 

4.2 

0.41 

4.3 

4.3 

0.41 

Table  9.  Volumetric  properties  of  numerical  model  sand/silty-clay  blends. 


Particle 
Fraction  % 
(coarse  to 
fine) 

15  psi 

50  psi 

Width  (in.) 

Height  (in.) 

Void  Ratio 

Width  (in.) 

Height  (in.) 

Void  Ratio 

0-100 

14.0 

14.2 

0.58 

13.5 

13.7 

0.44 

40-60 

3.3 

3.4 

0.58 

3.2 

3.3 

0.41 

48-52 

3.7 

3.7 

0.54 

3.5 

3.5 

0.41 

50-50 

3.8 

3.8 

0.49 

3.6 

3.6 

0.38 

52-48 

3.9 

3.9 

0.44 

3.7 

3.8 

0.36 

60-40 

4.2 

4.2 

0.40 

4.2 

4.2 

0.40 

The  bottom  wall  was  held  stationary,  simulating  the  pedestal  of  the  physi¬ 
cal  triaxial  cell.  The  top  wall  moved  toward  the  bottom  wall  at  a  fixed  rate 
under  velocity  control,  simulating  the  downward  movement  of  the  physi¬ 
cal  top  platen  under  the  force  of  the  axial  loading  piston.  To  maintain  the 
zero  volume  change  environment  of  undrained  triaxial  compression  test¬ 
ing,  for  each  increment  of  vertical  displacement,  the  lateral  confining  walls 
moved  outward  from  the  centroid  of  the  specimen  one-half  (1/2)  the  ver¬ 
tical  displacement  of  the  top  wall.  The  confining  walls  relaxed  outward 
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allowing  the  particulate  assembly  to  deform  both  vertically  and  radially, 
simulating  the  confinement  provided  by  the  neoprene  membrane  and 
applied  effective  confining  pressure  in  physical  triaxial  compression 
testing. 

The  simulated  testing  was  continued  for  a  fixed  number  of  computational 
increments  (200,000)  resulting  in  7.0  to  9.0  percent  shear  strain  in  the 
numerically  modeled  specimens.  Modeling  of  water  is  not  currently 
possible  using  the  DEM  model.  Therefore,  the  pore-water  pressure  (p)  of 
the  numerical  model  specimens  was  approximated  by  subtracting  the 
theoretical  total  mean  stress  from  the  effective  mean  stress  as  shown  in 
Equation  20. 


qi+2o3 

3 


(20) 


where: 

0!  =  principal  stress  measured  on  the  top  wall 

03  =  confining  pressure  calculated  as  the  average  stress  on  the  four 
lateral  confining  walls. 

Data  obtained  from  the  DEM  model  included  the  relative  position,  dis¬ 
placement,  force,  moment,  and  orientation  of  each  element  (spherical 
particles  and  confining  walls)  within  the  modeling  space  at  every  com¬ 
putational  increment. 

5.2  Testing  results 

Data  were  extracted  from  the  DEM  model  at  selected  intervals  to  include  1, 
10,  too,  1,000, 10,000,  50,000, 100,000, 150,000,  and  200,000  compu¬ 
tational  increments.  The  data  extracted  from  the  model  included  the  loca¬ 
tion  and  force  acting  on  each  of  the  six  confining  walls  and  each  of  the 
approximately  7,000  particles.  The  location  and  force  acting  on  the  con¬ 
fining  walls  were  used  to  calculate  the  shear  stress,  total  mean  stress, 
effective  mean  stress,  shear  strain,  and  pore-water  pressure.  The  confining 
walls  were  used  to  approximate  the  macroscale  properties  of  the  consoli¬ 
dated  specimen  that  would  traditionally  be  determined  in  physical  triaxial 
compression  testing.  The  location  and  forces  acting  on  the  individual  par- 
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tides  provided  data  not  available  in  standard  laboratory  testing  to  include 
the  distribution  of  stresses  within  discrete  particle-size  fractions  and  the 
evolution  of  coordination  number  during  shearing. 

5.2.1  Macroscale  mechanical  properties 

The  macroscale  mechanical  properties  of  each  specimen  determined  dur¬ 
ing  simulated  triaxial  compression  testing  included  shear  stress,  total 
mean  stress,  effective  mean  stress,  and  shear  strain,  calculated  using 
Equations  21  through  24,  respectively.  The  equation  for  determining  the 
pore-water  pressure  (Equation  20)  was  presented  previously. 

°s/!ear=0l-°3  (21) 


where: 


0 shear  =  shear  stress 

=  principal  stress  measured  on  the  top  wall 

03  =  confining  pressure  calculated  as  the  average  stress  on  the 
four  lateral  confining  walls. 


o 


o 


shear 


total  mean 


a. 


(22) 


where  Qlolai  mean  is  the  total  mean  stress. 


CT 


mean 


eg  +2a3 
3 


(23) 


where  0 effective  mean  is  the  effective  mean  stress. 

p  =15*p 

cshear  ±,yJ  caxial 


(24) 


where  &  shear  is  the  shear  strain  and  &  axial  is  the  axial  strain  measured  as 

the  change  in  distance  divided  by  the  initial  distance  between  the  top  and 
bottom  confining  walls. 
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5.2. 1.1  Model  calibration 

The  first  series  of  numerical  simulations  was  run  to  calibrate  the  DEM 
model.  The  characteristic  curves  developed  using  uniform  assemblies  of 
0.1875-in.  spheres  or  uniform  assemblies  of  sand/silty-clay,  as  presented 
in  Chapter  4,  were  used  to  calibrate  the  DEM  model.  For  each  material 
investigated  in  this  study  (stainless  steel,  polypropylene,  sand,  and  silty- 
clay)  the  material  properties  of  the  model  were  established  as  presented  in 
Table  6.  The  simulated  triaxial  compression  testing  was  run  for  200,000 
computational  increments  at  both  15-  and  50-psi  effective  confining  stress 
with  data  extracted  at  the  previously  presented  intervals.  The  extracted 
data  were  used  to  generate  stress-strain,  shear/principal  stress  (q-p),  and 
pore-water  pressure  plots.  The  data  from  the  simulated  testing  were  com¬ 
pared  to  the  characteristic  curves  developed  from  physical  laboratory 
testing  and  an  iterative  refinement  process  ensued  until  similar  mechani¬ 
cal  behavior  was  observed.  The  final  calibrated  material  properties  are 
presented  in  Table  10.  Comparisons  of  the  stress-strain  behavior  of  the 
physical  and  numerical  tests  are  presented  in  Figure  22  and  Figure  23  for 
the  ideal  spheres  and  naturally  occurring  materials,  respectively,  where 
the  physical  test  data  are  represented  with  a  solid  line  and  the  numerical 
test  data  are  represented  with  a  dashed  line. 


Table  10.  Calibrated  particle  properties  for  materials  simulated  in  DEM  model. 


Stainless  Steel 

Polypropylene 

Sand 

Silty-Clay 

Specific  gravity 

7.80 

0.90 

2.65 

2.70 

Normal  stiffness  (ksi) 

10,000 

1,500 

10,000 

1,200 

Shear  stiffness  (ksi) 

10,000 

1,500 

10,000 

1,200 

Rolling  resistance 

0.0 

0.0 

7.0 

0.0 

Torsion  resistance 

0.0 

0.0 

7.0 

0.0 

Translational  friction 

0.30 

0.35 

0.30 

0.50 

Rotational  friction 

0.30 

0.35 

0.30 

0.50 

5.2. 1.2  Model  validation 

Once  the  numerical  model  was  calibrated  using  the  physical  triaxial  com¬ 
pression  testing  data,  validation  of  the  model  was  pursued  by  examining 
bimodal  particle-size  distributions  similar  to  those  examined  in  the  labora¬ 
tory.  The  specimens  were  generated  and  tested  as  described  previously  in 
the  “Testing  program”  section.  Initially,  blends  of  40  percent  (by  volume) 
o.8-in.-diam  spheres  and  60  percent  o.2-in.-diam  spheres  were  examined 
at  both  15-psi  and  50-psi  effective  confining  pressure.  Macro-scale  behav- 
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ioral  plots  were  developed  including  stress-strain,  shear/mean  stress  (q-p 
plot),  and  pore-water  pressure  using  Equations  20  through  24.  This  process 
was  repeated  for  increasing  proportions  of  0.8-in.  diam  particles.  The 
designations  for  the  examined  particle-size  distributions  included  40-60, 
48-52,  50-50,  52-48,  and  60-40.  The  stress-strain  plots  for  the  numerically 
simulated  triaxial  compression  tests  run  in  this  investigation  at  15  and  50 
psi  are  presented  in  Figure  24  and  Figure  25  for  the  stainless  steel  blends, 
Figure  26  and  Figure  27  for  the  stainless  steel/polypropylene  blends,  and 
Figure  28  and  Figure  29  for  the  sand/silty-clay  blends. 


Figure  22.  Stress-strain  behavior  of  physical  (solid)  and 
model  (dashed)  uniform  spheres. 
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Figure  23.  Stress-strain  behavior  physical  (solid)  and  model  (dashed)  sand/silty-clay 
materials  from  Peters  and  Berney  (2010)  study. 
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Figure  24.  Stress-strain  behavior  of  simulated  stainless 
steel  spheres  at  15  psi. 
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Figure  25.  Stress-strain  behavior  of  simulated  stainless 
steel  spheres  at  50  psi. 
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Figure  26.  Stress-strain  behavior  of  simulated  polypropylene/ 
stainless  steel  spheres  at  15  psi. 
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Figure  27.  Stress-strain  behavior  of  simulated  polypropylene/ 
stainless  steel  spheres  at  50  psi 
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Figure  28.  Stress-strain  behavior  of  simulated  sand/ 
silty-clay  spheres  at  15  psi 
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Figure  29.  Stress-strain  behavior  of  simulated  sand/  silty-clay  spheres  at  50  psi. 
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In  examining  these  stress-strain  plots,  the  influence  of  confining  stress  is 
immediately  apparent.  At  low  confining  pressures  (15  psi),  the  differences 
in  mechanical  behavior  between  different  bimodal  gradations  is 
pronounced.  At  high  confining  pressures  (50  psi),  the  behaviors  become 
more  simular.  Trends  in  observed  mechanical  performance  appear  to  be 
more  greatly  influenced  by  confining  pressure  than  material  type. 
Regardless  of  composing  material,  at  15-psi  confining  pressure  the  40-60 
blend  displayed  the  weakest  shear  strength  and  at  50  psi,  the  60-40  blend 
appeared  the  weakest.  Clearly  defined  relationships  between  increased 
large  particle  proportion  and  shear  strength  at  both  high  and  low 
confining  pressure  levels  were  not  observed  for  any  of  the  bimodal 
materials  investigated.  When  compared  to  the  laboratory  triaxial  testing 
data  on  spherical  particles,  the  simulations  run  at  the  low  confining  stress 
(15  psi)  displayed  inverse  behavior.  However,  the  simulations  run  at  the 
high  confining  pressure  (50  psi)  displayed  similar  strength  trends  as  com¬ 
pared  to  the  laboratory  results.  Differences  in  expected  behavior  are  likely 
due  to  the  limited  number  of  particles  in  the  DEM,  which  allow  boundary 
conditions  to  influence  the  dilative  characteristics  within  the  simulation  at 
the  lower  confining  pressures. 
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5.2.2  Microscale  mechanical  properties 

The  primary  advantage  of  employing  numerical  testing  methods  is  the 
availability  of  particle  scale  data  not  easily  obtained  in  physical  testing. 
The  stresses  acting  on  and  connectivity  of  each  particle  within  the  com¬ 
pacted  assembly  are  data  available  only  through  employment  of  numerical 
testing  methods.  The  average  normal  and  shear  stress  acting  on  each  par¬ 
ticle  are  determined  by  summing  the  forces  acting  at  each  contact  that  a 
given  particle  participates  in  as  shown  in  Equation  25. 


(25) 


where: 

P 

Ojj  =  average  stress  tensor  component 
vp  =  volume  of  the  particle 

Nc  =  number  of  contacts  the  particle  participates  in 

Q 

Yj  =  distance  from  the  particle  centroid  to  the  point  of  contact 
fi  =  force  at  the  contact 


A  simplified  2-D  schematic  of  the  relationship  presented  in  Equation  25  is 
presented  in  Figure  30. 
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Figure  30.  Schematic  of  contact  force  components. 


Equation  25  is  used  to  describe  the  component  normal  and  shear  stresses 
that  define  the  stress  tensor  for  each  particle.  The  component  stresses 
were  used  in  this  study  to  calculate  the  effective  mean  and  shear  stresses 
for  each  particle  in  the  assembly  as  shown  in  Equations  26  and  27. 
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effective  mean 


°xx  +  °yy  +  °zz 

3 


(26) 


where 
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effective  mean 
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XX 


a 


yy 


a 


zz 


effective  mean  stress 
normal  stress  in  the  x-direction 
normal  stress  in  the  y-direction 
normal  stress  in  the  z-direction 


a 


shear 


where: 


a 


shear 


(27) 


=  shear  stress 
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°xy  =  x~y  shear  stress  component 
G xz  =  x-z  shear  stress  component 
°yz  =  y~z  shear  stress  component 

5.2.2. 1  Distribution  of  stresses 

The  distribution  of  stress  amongst  the  particle-size  fractions  of  a  bimodal 
blend  provides  insight  into  the  proportion  that  is  dominating  engineering 
behavior.  The  normal  and  shear  stresses  were  averaged  for  each  particle 
size  for  the  bimodal  blends  investigated  in  this  study.  The  comparisons  of 
average  stress  within  each  particle-size  fraction  for  the  stainless  steel, 
polypropylene/stainless  steel,  and  sand/ silty-clay  numerical  specimen 
after  consolidation  at  both  15-  and  50-psi  effective  confining  pressure  are 
presented  in  Figure  31  through  Figure  33,  respectively. 


Figure  31.  Shear  stress  distribution  for  polypropylene/stainless  steel 
specimens  at  15  psi  (left)  and  50  psi  (right)  confining  pressure. 


Figure  32.  Shear  stress  distribution  for  stainless  steel  specimens  at 
15  psi  (left)  and  50  psi  (right)  confining  pressure. 
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Figure  33.  Shear  stress  distribution  for  sand/silty  clay  specimens  at 
15  psi  (left)  and  50  psi  (right)  confining  pressure. 


Clear  trends  when  examining  the  distribution  of  stress  amongst  the  dis¬ 
crete  particle-size  fractions  are  not  observed  following  isotropic  consolida¬ 
tion  (Figures  32  and  33).  A  state  of  stress  equilibrium  is  approached  at 
both  15  and  50  psi  for  the  stainless  steel  assemblies  blended  at  or  near  50- 
50  and  the  polypropylene/stainless  steel  assemblies  blended  at  48-52.  No 
stress  equilibrium  state  is  observed  for  the  bimodal  sand/silty  clay  blends 
investigated  in  this  study.  During  shear,  regardless  of  particle-size 
distribution  or  composing  material  type,  the  majority  of  the  stress  is 
transmitted  through  the  small  size  fraction. 

5. 2.2.2  Coordination  number 

Coordination  number  is  simply  a  measure  of  the  number  of  contacts  a 
given  particle  participates  in.  The  average  coordination  number  for  a  spe¬ 
cimen  is  the  total  number  of  contacts  in  the  assembly  divided  by  the  num¬ 
ber  of  particles.  The  average  coordination  number  provides  an  indication 
of  the  connectivity  of  the  assembly  and  also  provides  an  indication  of  the 
development  of  force  chains.  An  increase  in  average  coordination  number 
is  commonly  observed  during  consolidation.  As  the  particles  are  forced 
closer  together  the  number  of  contacts  increases.  During  dilation,  the 
specimen  is  increasing  in  volume,  reducing  the  density  of  particles  and 
subsequently,  the  coordination  number.  As  force  chains  develop,  the  coor¬ 
dination  number  of  the  participating  particles  increases,  while  the  coordi¬ 
nation  number  of  particles  not  participating  in  the  force  chain  decreases. 
The  average  total  coordination  number,  average  large  proportion  coordi¬ 
nation  number,  and  average  small  proportion  coordination  number  fol¬ 
lowing  isotropic  consolidation  for  the  specimens  examined  in  this  study  at 
both  15-  and  50-psi  effective  confining  pressure  are  presented  in  Tables  11 
through  13. 
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Table  11.  Coordination  number  for  stainless  steel  specimens. 


Particle 
Fraction  % 
(coarse  to 
fine) 

15  psi 

50  psi 

Total 

0.80-in. 

Particles 

0.20-in. 

Particles 

Total 

0.80-in. 

Particles 

0.20-in. 

Particles 

0-100 

6.2 

- 

6.2 

6.6 

- 

6.6 

40-60 

6.6 

30.1 

6.4 

8.5 

37.2 

8.2 

48-52 

7.3 

36.8 

6.9 

8.5 

43.1 

8.0 

50-50 

7.6 

39.4 

7.1 

8.6 

44.3 

8.0 

52-48 

6.7 

33.6 

6.2 

8.8 

45.4 

8.2 

60-40 

8.0 

39.9 

7.2 

8.1 

40.3 

7.3 

Table  12.  Coordination  number  for  polypropylene/ 
stainless  steel  specimens. 


Particle 
Fraction  % 
(coarse  to 
fine) 

15  psi 

50  psi 

Total 

0.80-in. 

Particles 

0.20-in. 

Particles 

Total 

0.80-in. 

Particles 

0.20-in. 

Particles 

0-100 

6.2 

- 

6.2 

6.6 

- 

6.6 

40-60 

6.8 

30.3 

6.5 

8.5 

37.5 

8.2 

48-52 

7.7 

38.6 

7.3 

8.8 

44.1 

8.3 

50-50 

7.8 

40.3 

7.3 

8.5 

43.9 

8.0 

52-48 

8.8 

45.4 

8.2 

8.9 

45.6 

8.3 

60-40 

8.2 

41.1 

7.4 

8.2 

41.2 

7.4 

Table  13.  Coordination  number  for  sand/silty  clay  specimens. 


Particle 
Fraction  % 
(coarse  to 
fine) 

15  psi 

50  psi 

Total 

0.80-in. 

Particles 

0.20-in. 

Particles 

Total 

0.80-in. 

Particles 

0.20-in. 

Particles 

0-100 

7.1 

- 

7.1 

8.2 

-- 

8.2 

40-60 

8.1 

36.5 

7.8 

9.2 

40.0 

8.9 

48-52 

8.5 

43.8 

8.1 

9.6 

48.2 

9.1 

50-50 

8.4 

43.4 

7.9 

9.5 

46.5 

8.9 

52-48 

8.7 

45.1 

8.1 

9.4 

48.0 

8.7 

60-40 

8.6 

42.9 

7.8 

8.7 

43.1 

7.8 

As  a  general  trend,  total,  large  proportion,  and  small  proportion  coordina¬ 
tion  number  increases  with  increasing  large  proportion  content.  Regard¬ 
less  of  material  and  confining  pressure  level,  the  52-48  blends  appear  to 
have  the  largest  total,  large  proportion,  and  small  proportion  coordination 
numbers.  During  shear  the  coordination  number  for  all  of  the  investigated 
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assemblies  decreases,  with  the  greatest  decrease  occurring  in  the  60-40 
blends. 
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6  Analysis  and  Discussion 

6.1  Macroscopic  behavior 

Similar  to  the  observations  of  O’Sullivan  et.  al  (2006),  the  DEM  model  was 
slightly  stiffer  and  reached  peak  strength  at  a  lower  level  of  shear  strain  as 
compared  to  laboratory  data.  However,  the  DEM  model  successfully 
captured  volumetric  strain  response.  O’Sullivan  et.  al  (2004)  observed  that 
as  distribution  of  particle  sizes  increases,  peak  strength  decreases  and 
peak  strain  at  failure  increases.  A  definite  variation  in  the  response  curve 
was  observed  to  be  associated  with  the  standard  deviation  of  the  particle- 
size  distribution  (O’Sullivan  et.  al  2006).  For  the  bimodal  distributions 
investigated  in  this  study  (Figure  22  through  Figure  29),  greater  shear 
strength  and  decreased  shear  strain  were  observed  for  the  blends  as 
compared  to  the  uniform  assemblies,  with  the  exception  of  the  sand/silty 
clay  blends  at  50-psi  confining  pressure  that  fell  within  the  large  and  small 
particle  uniform  assemblies  (Figure  29).  A  definite  variation  in  response 
was  observed  for  the  specimens  tested  at  the  low  confining  pressure 
(15  psi),  but  the  response  was  more  greatly  influenced  by  the  coordination 
number  than  the  distribution  of  particle  sizes. 

Cui  et.  al  (2007)  observed  that  polydisperse  specimens  exhibited  greater 
radial  dilation  as  compared  to  monodisperse  specimens.  The  same 
observation  was  made  in  this  study.  The  uniform  materials  investigated  in 
this  study  including  the  stainless  steel,  polypropylene,  sand,  and  silty-clay 
showed  reduced  dilative  tendency  as  compared  to  the  bimodal  blends.  The 
less  stiff  materials,  polypropylene  and  silty-clay,  exhibited  purely 
contractive  behavior,  while  the  stiffer  materials,  stainless  steel  and  sand, 
exhibited  a  short  period  of  contraction  followed  by  dilation.  All  of  the 
investigated  bimodal  blends  exhibited  purely  dilative  behavior. 

Cavarretta  et.  al  (2009)  observed  that  overall  response  is  dominated  by  the 
particles  participating  in  force  chains.  An  in-depth  investigation  into  the 
formation  and  evolution  of  force  chains  was  not  conducted  as  part  of  this 
study.  However,  the  theory  of  percolation  provides  that  once  the 
percolation  threshold  is  reached,  sufficient  coarse  particles  exist  for  the 
formation  of  force  chains.  Nearing  the  percolation  threshold  is  indicated 
by  erratic  behavior  and  a  reduction  in  coordination  number.  In  this 
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investigation,  a  clear  transition  in  dominant  particle  behavior  was  not 
observed.  Increases  in  coarse  material  content  resulted  in  gradual  changes 
in  mechanical  behavior.  However,  when  examining  the  sand/silty-clay 
blends  tested  at  the  low  confining  pressure  (15  psi),  it  can  be  observed  that 
as  coarse  particle  content  increases  the  stress-strain  curve  shifts  between 
the  uniform  fine  and  uniform  coarse  assemblies  (Figure  28). 

The  macro-scale  response  is  a  function  of  particle  shape,  particle  size,  and 
void  distributions  in  addition  to  other  internal  variables  (Thornton  2000). 
Macro-scale  response  includes  deformability,  strength,  dilatancy,  and 
strain  localization  (Belheine  et  al.  2009).  The  influences  of  particle-size 
distribution  and  particle  stiffness  were  examined  in  this  study.  Distinct 
changes  in  response  resulting  from  changes  in  particle-size  distribution 
were  not  observed.  However,  coordination  number  was  found  to  highly 
influence  material  behavior.  Additionally,  the  properties  of  the  material 
were  found  to  have  a  less  significant  influence.  The  blends  composed  of 
stainless  steel  only  and  the  blends  composed  of  less  stiff  polypropylene 
and  stainless  steel  were  observed  to  have  similar  mechanical  response. 
Common  coordination  numbers  were  measured  for  these  blends. 

6.2  Microscopic  behavior 

Stresses  calculated  from  the  boundary  forces  were  smaller  than  the 
stresses  measured  internally  due  to  non-uniformities  within  the  specimen 
(Cui  et.  al  2007).  The  bimodal  blends  examined  in  this  study  exhibited 
contrary  behavior.  The  stresses  calculated  from  the  bounding  walls  were 
found  to  be  significantly  larger  than  the  stresses  measured  within  the 
particles.  This  is  hypothesized  to  be  due  to  the  use  of  rigid  confining  walls 
versus  the  flexible  membrane  employed  by  Cui  et.  al  (2007).  In  the  tightly 
compacted  specimen,  the  resultant  forces  measured  at  the  rigid  confining 
walls  can  be  influenced  by  a  single  highly  stressed  particle.  Thus,  the 
boundary  stress  measured  at  the  wall  can  be  significantly  larger  than  the 
average  stress  of  the  particulate  assembly. 

Micro-scale  properties  influencing  macro-scale  response  include  normal, 
tangential,  and  rolling  stiffness,  local  friction,  and  non-dimensional  plastic 
coefficient  at  the  contacts  (Belheine  et  al.  2009).  Particle  stiffness,  assem¬ 
bly  fabric,  and  shear  history  determine  the  distribution  of  forces  within  a 
particulate  system  (Antony  2000).  The  influence  of  particle  properties  to 
include  specific  gravity,  contact  stiffness,  surface  friction,  and  rolling 
resistance  were  examined  in  this  study.  Values  of  specific  gravity  ranged 
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from  0.9  to  7.8,  values  of  contact  stiffness  ranged  from  10,000  to  1,200, 
values  of  surface  friction  ranged  from  0.3  to  0.5,  and  a  rolling  resistance 
value  of  7.0  was  applied  to  the  simulated  sand  particles  to  account  for  the 
use  of  spherical  particles  to  model  non-spherical  materials.  All  of  the 
applied  properties  greatly  influenced  the  behavior  of  the  modeled 
materials  and  enabled  reasonable  replication  of  laboratory  testing  of 
physical  specimens. 

Forces  are  transmitted  in  well-graded  specimens  through  systems  of  force 
chains  versus  well  distributed  in  uniformly  graded  specimens  (Evans  et  al. 
2009).  Antony  (2000)  observed  that  normal  contact  forces  can  be  subdi¬ 
vided  into  those  carrying  less  than  the  average  and  those  carrying  more. 
While  these  findings  are  well  accepted  for  physical  specimens,  these 
behaviors  were  not  observed  in  this  study.  The  selection  of  bimodal 
particle  sizes,  0.80-  and  o.20-in.-diameter  particles,  and  the  specimen  size 
did  not  allow  for  the  formation  of  well-established  force  chains.  The  1:4 
large  particle  to  small  particle  diameter  size  ratio  was  not  sufficient  to 
allow  the  small  proportion  particles  to  fit  into  the  void  between  contacting 
large  particles.  This  prevented  the  short-circuiting  of  the  small  proportion, 
requiring  the  entire  particle  assembly  to  share  the  applied  stresses. 
Additionally,  the  large  particle  diameter  to  specimen  size  ratio  of  1:5  was 
not  sufficient  to  allow  for  the  formation  of  force  chains. 

6.2.1  Distribution  of  stresses 

The  distribution  of  shear  stresses  presented  in  Figure  31  through  Figure  33 
show  the  average  shear  stress  within  each  particle-size  fraction.  However, 
to  understand  the  contribution  of  each  particle-size  fraction  to  the 
strength  capability  of  the  entire  specimen,  the  volume  of  each  particle 
must  also  be  considered.  To  capture  the  contribution  of  each  particle 
fraction  to  the  entire  assembly  strength,  the  stress  porosity  was 
considered.  In  this  technique  the  stress  within  each  particle  is  weighted 
based  upon  its  volume  according  to  Equation  28  and  averaged  across  each 
particle-size  fraction.  Figure  34  through  Figure  36  show  the  stress  porosity 
for  the  stainless  steel,  polypropylene/stainless  steel,  and  sand/silty-clay 
bimodal  assemblies  respectively.  A  clear  transition  in  material  behavior  is 
not  apparent  for  the  assemblies  containing  the  spherical  particles. 
However,  the  assemblies  of  sand/silty-clay  show  increasing  stress  porosity 
for  the  large  material  as  the  proportion  is  increased. 
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Figure  34.  Stress  porosity  of  stainless  steel  specimens  at  15  psi  (left) 
and  50  psi  (right)  confining  pressure. 


Figure  35.  Stress  porosity  of  polypropylene/stainless  steel  specimens  at 
15  psi  (left)  and  50  psi  (right)  confining  pressure. 


Figure  36.  Stress  porosity  of  sand/silty  clay  specimens  at  15  psi  (left) 
and  50  psi  (right)  confining  pressure. 
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6.2.2  Coordination  number 

Coordination  number  describes  the  average  number  of  contacts  for  the 
particles  in  the  assembly  (Bathurst  and  Rothenburg  1988).  Particulate 
materials  transmit  loads  through  shared  contacts.  Coordination  number 
strongly  influences  the  behavior  of  particulate  materials  because  increased 
coordination  number  results  in  more  uniform  distribution  of  loads.  O’Sul¬ 
livan  et.  al  (2006)  observed  decreased  strength  accompanied  a  reduction 
in  average  coordination  number.  The  critical  coordination  number  reflects 
a  minimum  requirement  for  physical  stability  (Thornton  2000).  A 
coordination  number  of  approximately  6,  for  dense  assemblies,  or  4,  for 
loose  assemblies,  indicates  the  required  contacts  for  stability  have  been 
achieved  and  the  system  effectively  locks  up.  When  the  critical 
coordination  number  is  reached,  significant  volume  change  is  observed. 
Coordination  number  greatly  influences  shear  displacement  with 
increasing  number  of  contacts  limiting  particle  rotation.  Assemblies  with 
coordination  numbers  less  than  3  are  unstable.  As  coordination  number 
reaches  its  maximum  value  of  6,  particles  are  restrained  to  a  point  that 
rotations  practically  disappear  (Bathurst  and  Rothenburg  1988). 

The  critical  coordination  number  described  by  Bathurst  and  Rothenburg 
(1988)  describes  uniform  assemblies  of  rigid  spheres.  The  coordination 
number  for  the  uniform  assemblies  of  spheres  examined  in  this  study 
exceeded  the  proposed  maximum  value  of  6.  Cui  et.  al  (2007)  indicated 
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that  the  initial  coordination  number  for  polydisperse  specimens  is  lower 
than  monodisperse  specimens  despite  lower  void  ratio.  This  relationship 
was  not  observed  in  this  study.  The  void  ratio  for  the  bimodal  numerical 
specimens  generated  as  part  of  this  study  were  lower  that  uniform 
assemblies,  but  also  exhibited  greater  initial  coordination  number. 

During  triaxial  compression  testing,  the  coordination  number  decreased 
slightly  for  the  bimodal  distributions,  but  more  significantly  for  the  uni¬ 
form  assemblies.  Although  the  specimen  coordination  number  stabilizes 
after  3  percent  strain,  the  coordination  number  within  the  specimen 
continually  changes  with  some  regions  increasing  and  others  decreasing 
(Cui  et.  al  2007).  After  an  initial  rearrangement,  the  particle  matrix  in 
granular  systems  remains  constant.  Shearing,  or  change  in  stress  ratio, 
results  in  reduced  coordination  number  (Sitharam  et.  al  2002).  The 
change  in  coordination  number  of  the  stainless  steel, 
polypropylene/stainless  steel,  and  sand/silty-clay  bimodal  specimens 
between  o  and  5  percent  shear  strain  is  presented  in  Tables  14  through  16, 
respectively. 

During  shearing,  particles  contributing  in  force  chains  undergo  continuous 
cycles  of  compression  and  relaxation  without  a  change  in  coordination 
number  (Antony  2000).  Although  small  changes  in  coordination  number 
for  the  assembly  are  measured,  particles  participating  in  force  chains  may 
have  increasing  coordination  numbers  that  are  offset  by  the  reduced 
coordination  number  of  the  non-contributing  particles  in  the  assembly. 
O’Sullivan  et.  al  (2004)  observed  that  traditional  measures  of  fabric  such 
as  coordination  number  are  not  effective  for  distinguishing  differences  in 
mechanical  performance.  While  coordination  number  may  not  be  effective 
for  determining  the  formation  and  evolution  of  force  chains,  results  of  this 
study  indicate  that  total  coordination  number  correlates  well  with 
observed  mechanical  performance.  For  similar  composing  materials, 
increased  values  of  coordination  number  corresponded  with  increased 
shear  strength.  The  shear  stress  and  coordination  number  at  5  percent 
shear  strain  were  plotted  for  the  stainless  steel,  polypropylene/stainless 
steel,  and  sand/silty-clay  bimodal  blends  investigated  in  this  study  in 
Figure  37  and  Figure  38.  A  clear  linear  trend  between  total  coordination 
number  and  shear  stress  is  apparent. 
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Table  14.  Change  in  coordination  number  during  shear 
of  stainless  steel  specimens. 


Particle 
Fraction  % 
(coarse  to 
fine) 

15  psi 

50  psi 

CN 

0%  Eshear 

CN 

5%  Eshear 

AON 

CN 

0%  Eshear 

CN 

5%  Eshear 

ACN 

0-100 

6.2 

5.7 

-8.9% 

6.6 

5.9 

-10.6% 

40-60 

6.6 

6.5 

-0.9% 

8.5 

8.5 

-0.6% 

48-52 

7.3 

7.2 

-0.4% 

8.5 

8.5 

-0.3% 

50-50 

7.6 

7.6 

-0.5% 

8.6 

8.5 

-1.1% 

52-48 

6.7 

6.6 

-1.1% 

8.8 

8.7 

-1.0% 

60-40 

8.0 

7.7 

-3.6% 

8.1 

7.8 

-3.6% 

Table  15.  Change  in  coordination  number  during  shear  of 
polypropylene/stainless  steel  specimens. 


Particle 
Fraction  % 
(coarse  to 
fine) 

15  psi 

50  psi 

CN 

0%  Eshear 

CN 

5%  Eshear 

ACN 

CN 

0%  Eshear 

CN 

5%  Eshear 

ACN 

0-100 

6.2 

5.7 

-8.9% 

6.6 

5.9 

-10.6% 

40-60 

6.8 

6.6 

-1.8% 

8.5 

8.5 

-0.2% 

48-52 

7.7 

7.7 

-0.1% 

8.8 

8.7 

-0.8% 

50-50 

7.8 

7.7 

-0.5% 

8.5 

8.5 

-0.5% 

52-48 

8.8 

8.7 

-1.3% 

8.9 

8.8 

-1.2% 

60-40 

8.2 

7.9 

-3.5% 

8.2 

7.9 

-3.4% 

Table  16.  Change  in  coordination  number  during  shear 
of  sand/silty  clay  specimens. 


Particle 
Fraction  % 
(coarse  to 
fine) 

15  psi 

50  psi 

CN 

0%  Eshear 

CN 

5%  Eshear 

ACN 

CN 

0%  Eshear 

CN 

5%  Eshear 

ACN 

0-100 

7.1 

6.2 

-13.9% 

8.2 

8.2 

0.0% 

40-60 

8.1 

8.1 

-0.2% 

9.2 

9.1 

-0.9% 

48-52 

8.5 

8.5 

-0.2% 

9.6 

9.5 

-0.8% 

50-50 

8.4 

8.4 

-0.4% 

9.5 

9.4 

-1.1% 

52-48 

8.7 

8.6 

-0.7% 

9.4 

9.3 

-1.0% 

60-40 

8.6 

8.3 

-3.0% 

8.7 

8.4 

-2.9% 
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Figure  37.  Influence  of  coordination  number  on  shear  strength  of 
specimens  sheared  at  15-psi  confining  pressure. 
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Figure  38.  Influence  of  coordination  number  on  shear  strength  of 
specimens  sheared  at  50-psi  confining  pressure. 
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7  Summary  and  Conclusions 

The  influence  of  particle-size  distribution  on  the  stress  transfer  mecha¬ 
nisms  occurring  in  bimodal  particulate  materials  was  investigated  in  this 
study.  A  serial  version  of  the  3-D  DEM  model  developed  by  Horner  et.al 
(2001)  was  used  to  simulate  the  consolidated-undrained  triaxial 
compression  testing  of  cubical  specimens.  The  DEM  model  was  calibrated 
using  laboratory  tests  on  assemblies  of  uniform  spheres  composed  of 
stainless  steel,  polypropylene,  sand,  and  silty-clay  materials.  Calibrated 
material  properties  were  developed  by  comparing  the  stress-strain 
behavior  of  simulated  tests  to  the  results  obtained  in  the  laboratory.  Once 
calibrated,  the  model  was  used  to  simulate  triaxial  testing  on  bimodal 
specimens  composed  of  0.80-  and  o. 20-in. -diam  spheres.  The  stress- 
strain  response  of  specimens  subjected  to  simulated  triaxial  compression 
testing  compared  favorably  to  the  relationships  observed  in  laboratory 
testing. 

The  primary  advantage  to  employing  numerical  testing  techniques  as 
compared  to  traditional  laboratory  testing  is  the  availability  of  micro-  or 
particle-scale  data.  Data  obtained  from  the  DEM  in  this  study  included 
individual  particle  and  particle-size  fraction  average  stresses  and 
coordination  number.  Distinct  changes  in  behavior  resulting  from  changes 
in  particle-size  fractions  were  not  observed.  However,  a  relationship 
between  coordination  number  and  specimen  stiffness  was  observed. 
Coordination  number  describes  the  average  number  of  contacts  for  each 
particle  within  a  compacted  assembly.  Coordination  number  is  influenced 
by  a  number  of  factors  to  include  particle  shape,  surface  friction,  size  ratio, 
and  assembly  density.  The  density  is  strongly  influenced  by  the  proportion 
of  particle  sizes  composing  the  assembly  and  therefore,  particle-size 
fraction  has  an  indirect  effect  on  the  shear  strength  of  particulate 
materials. 

The  central  theme  of  this  work  was  related  to  percolation  theory  as 
described  in  Peters  and  Berney  (2010).  Percolation  theory  describes  the 
critical  connectivity  within  a  network.  The  point  at  which  the  critical  con¬ 
nectivity  is  reached  is  the  percolation  threshold.  With  respect  to  bimodal 
particulate  materials,  the  exhibited  behavior  is  controlled  by  the  dominant 
size  fraction.  As  the  subordinate  fraction  is  increased,  the  percolation 
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threshold  is  reached,  and  the  behavior  transitions  from  being  dominated 
by  one  fraction  to  the  other.  Peters  and  Berney  (2010)  observed  erratic 
and  even  unstable  behavior  in  blends  of  sand  and  silty-clay  at  the  percola¬ 
tion  threshold.  Once  the  percolation  threshold  is  reached  in  a  fine/coarse 
bimodal  system,  sufficient  connectivity  exists  for  the  formation  of  force 
chains.  The  development  of  force  chains  results  in  stress  localization  and 
short-circuiting  of  the  fine  proportion. 

The  percolation  threshold  was  investigated  in  this  study  by  examining 
bimodal  distributions  that  were  dominated  by  one  fraction  or  were 
approaching  the  percolation  threshold.  The  investigated  blends  included 
0-100,  40-60,  48-52,  50-50,  52-48,  and  60-40  (coarse-fine)  content  by 
volume  percentage.  Additionally,  the  influence  of  material  type  on  the 
critical  particle-size  proportions  was  investigated  by  examining  blends  of 
stainless  steel  only,  less  stiff  large  fraction  composed  of  polypropylene 
blended  with  small  stainless  steel,  and  less  stiff  small  fraction  in  the  blends 
of  sand  and  silty-clay  as  presented  in  Peters  and  Berney  (2010).  Distinct 
changes  in  behavior  resulting  from  changes  in  particle-size  fraction  were 
not  observed  in  this  study  with  the  exception  of  the  sand/silty-clay  blends 
tested  at  a  confining  pressure  of  15  psi.  The  stress-strain  behavior  of  these 
blends  transitioned  from  resembling  a  uniform  assembly  of  fine  particles 
to  a  uniform  assembly  of  coarse  particles  with  increasing  coarse  material 
proportion.  Although  no  definite  transitional  proportion  was  identified, 
the  theory  of  percolation  was  validated  for  the  particulate  materials 
observed  in  this  study.  The  objective  of  this  investigation  was  to  examine 
the  influence  of  particle-size  distribution  on  the  stress  transfer 
mechanisms  of  bimodal  particulate  assemblies.  Both  physical  and 
numerical  methods  were  employed.  The  physical  specimens  showed 
definite  changes  in  stress-strain  behavior  resulting  from  changes  in 
particle-size  distribution.  However,  the  expected  increases  and  decreases 
in  strength  based  on  large  particle  proportion  was  not  observed  principally 
due  to  the  limiting  number  of  particles  used  within  the  DEM  simulation 
which  allowed  boundary  conditions  to  influence  the  overall  constitutive 
response.  The  stress-strain  behavior  of  the  numerical  specimens  differed 
both  from  expectation  and  from  the  observed  laboratory  performance. 
However,  a  relationship  between  coordination  number  and  assembly 
stiffness  was  observed.  Particulate  materials  transmit  loads  at  the  shared 
contact  points.  For  the  bimodal  assemblies  investigated  in  this  study,  an 
increase  in  contacts  resulted  in  greater  sharing  of  loads  as  indicated  by 
greater  strength  and  stiffness,  resulting  in  more  efficient  stress  transfer. 
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Factors  contributing  to  increased  contacts  or  coordination  number  can 

also  contribute  to  increased  stiffness  and  strength. 

Key  findings  from  this  investigation  include: 

•  DEM  is  an  effective  tool  for  investigating  the  micromechanics  driving 
macro-scale  performance  in  particulate  materials. 

•  Careful  calibration  is  required  to  model  materials  of  interest  using 
DEM  models. 

•  Differences  in  observed  behavior  between  physical  and  numerical 
experimentation  methods  are  owing  to  the  formulation  of  the 
numerical  model.  While  not  all  issues  can  be  addressed,  careful 
modeling  of  the  physical  system  results  in  good  agreement  between 
physical  and  numerical  methods. 

•  The  existence  of  a  distinct  percolation  threshold  was  not  observed  in 
this  investigation. 

•  The  behavior  of  both  physical  and  numerical  specimens  with  changing 
particle-size  fractions  differed  from  expectations. 

•  While  no  direct  relationship  between  particle-size  distribution  and 
macro-scale  behavior  were  observed,  a  distinct  relationship  between 
coordination  number  and  stress-strain  behavior  was  documented. 

•  The  domain  size  of  the  DEM  of  7,000  particles  was  likely  not  sufficient 
to  elicit  the  percolation  response  in  the  particle  systems.  A  larger 
domain  size  would  reduce  boundary  influences  and  allow  for  a  clearer 
evaluation  of  particle  response. 

Recommendations  resulting  from  this  work  include: 

•  Further  investigation  is  required  to  quantify  the  influence  of  modeled 
features  to  include  rigid  walls  and  cubical  specimens. 

•  Further  development  of  the  Horner,  Peters,  and  Carrillo  (2010)  model 
is  required  to  include  the  use  of  flexible  membrane  and  pore  water 
fluid. 

•  Additional  testing  is  required  on  non-spherical  particles  to  examine  the 
validity  of  the  coordination  number  and  stress-strain  behavior  rela¬ 
tionship  in  naturally  occurring  particles. 

•  Use  of  a  larger  particle  domain  size  to  enable  percolation  responses  to 
develop.  This  will  require  greater  computer  processing  power  and 
parallelization  of  the  DEM  code  to  optimize  run  times. 
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